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.

Improve this question
$\endgroup$
4
  • 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
    Commented Nov 15, 2019 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
    Commented Nov 15, 2019 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
    Commented Nov 15, 2019 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.'s missing motivation ♦
    Commented Nov 15, 2019 at 22:36

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
1
  • $\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
    Commented Nov 16, 2019 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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