0

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.

davideor
  • 1
  • 2
  • 1
    To get 500 random numbers from your distribution use SeedRandom[1234]; 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 https://mathematica.stackexchange.com/questions/13981/how-to-generate-a-randomvariate-of-a-custom-distribution – 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.'s missing motivation Nov 15 '19 at 22:36

1 Answers1

1

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.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
  • 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