0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ To get 500 random numbers from your distribution use SeedRandom[1234]; data = RandomVariate[dist, 500]; The SeedRandom is needed to give you reproducible results. $\endgroup$ – Bob Hanlon Nov 15 '19 at 20:04
  • $\begingroup$ 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. $\endgroup$ – davideor Nov 15 '19 at 20:25
  • $\begingroup$ 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/… $\endgroup$ – davideor Nov 15 '19 at 20:47
  • $\begingroup$ 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. $\endgroup$ – J. M. is in limbo Nov 15 '19 at 22:36
1
$\begingroup$

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[42]
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.

$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – davideor Nov 16 '19 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.