# Generating random numbers from specified distribution under a constraint

I need to generate 500 numbers from a distribution specified. I used the well-known CDF method, but since I'm new to Mathematica, I don't know how to set a condition on the sum of these variables.

dist :=
ProbabilityDistribution[(2/Sqrt[Pi]* Exp[-1/x])/x^(5/2), {x, 0, ∞}];
g[x_] = CDF[dist, x];
h =
Table[{g[x], x}, {x, 0, 5500}] // N // Interpolation[#, InterpolationOrder -> 5] &;

d = 1;
While[d < 5990 || d > 6010,
h /@ (RandomReal[{0, 1}, {500}]);
d = Total[h]]


Where, with Total, I'm trying to obtain the sum of the interpolated array but for some reason it gives me this error:

Lists of unequal length in InterpolatingFunction...

Update:

d = 1; While[d < 5990 || d > 6010,
a = h /@ (RandomReal[{0, 1}, {500}]);
d = Total[a]]


This works. Tahnks anyway for the help.

• To get 500 random numbers from your distribution use SeedRandom; data = RandomVariate[dist, 500]; The SeedRandom is needed to give you reproducible results. – Bob Hanlon Nov 15 '19 at 20:04
• Probably I didn’t explain myself properly: I need 500 values from my distribution conditioned by the value of their sum. For this reason I need a cycle until I have the set I’m looking for. – davideor Nov 15 '19 at 20:25
• Moreover the solution you suggested gives quite different results compared to the one here suggested, which is the one I used mathematica.stackexchange.com/questions/13981/… – davideor Nov 15 '19 at 20:47
• FWIW, your CDF has a closed form: 2/Sqrt[π t] Exp[-1/t] + Erfc[1/Sqrt[t]]. You can perhaps try using it in FindRoot[] to implement the inverse CDF. Alternatively, you could look into using NDSolve[], as was done here. – J. M. is in limbo Nov 15 '19 at 22:36

Your probability distribution has a mean of 2. That means the expected value of 500 random variates will be 1000. Therefore, the following will work:

SeedRandom
sample =
With[{sum = 1000, d = 10},
Module[{testSum = 1, sample},
While[testSum < sum - d || testSum > sum + d,
sample =
Quiet @
RandomVariate[
ProbabilityDistribution[
(2/Sqrt[Pi] Exp[-1/x])/x^(5/2), {x, 0, ∞}],
500];
testSum = Total[sample]];
sample]];

Short[sample, 4]


{1.07783, 1.24783, 0.550867, 20.2222, 0.185783 <<491>>, 0.923722, 1.62345, 0.47509, 0.727451}

Total[sample]


993.243

However, if you were to raise the value of sum to 6000, your chance drawing a sample meeting your requirement would less than the chance of a snowball successfully rolling through hell.

• I suggest you trying the code in the update, it's in fact an interesting phenomenon and it's named "condensation". (The simulation lasts a couple of seconds) – davideor Nov 16 '19 at 10:37