3
$\begingroup$

I need to generate a random number (real) and then test if this number works on my pdf. if it works i append it to a list. if it doesn't, i reject it.

Until now i have done this:

f[x_] := x/2; "with 0<=x<=2"

this is the function i'm using

testexc = {};
xc = {};

These are the lists where i want to append numbers

While[Length[testexc] <= 10000, 
  AppendTo[testexc, RandomReal[{0, 1}]]];

This is how i generate 10000 Random Numbers within the specified range

What i want now is: I need 10000 random numbers that works with my function. And i couldnt find out any way of doing it.

$\endgroup$
  • $\begingroup$ There is standard answer 2635 to such questions. $\endgroup$ – garej Jul 10 at 8:33
2
$\begingroup$

Inverse CDF method

Here, try this; it should be faster:

pdf = x \[Function] x/2;
cdf = x \[Function] Evaluate[Integrate[pdf[t], {t, 0, x}]]
cdfinv = y \[Function] 2 Sqrt[y]
rand = cdfinv[RandomReal[{0, 1}, {1000000}]]; // RepeatedTiming // First

0.009

One million random number is a percent of a second.

Plotting a histogram to check the distribution is correct:

Histogram[rand]

enter image description here

Actually, there is also a built-in method for this. It goes like this:

distro = ProbabilityDistribution[x/2, {x, 0, 2}];
rand = RandomVariate[distro, 1000000]; // RepeatedTiming // First

0.08

For some reason, it is significantly slower...

Acceptance/rejection method

A listable approach

If you insist on using the "acceptance/rejection method" (better know as Monte Carlo method, you can do this:

n = 2000000;
First@RepeatedTiming[

  x = RandomReal[{0, 2}, n];
  y = RandomReal[{0, 1}, n];
  rand = Pick[x, UnitStep[Subtract[pdf[x], y]], 1];

  ]

0.062

This generates about a million random numbers with 0.062 seconds. I would strongly discourage methods that use Append repeatedly, because the will quadratic complexity and are very memory bound (each time Append is called, you have to copy the full array).

An approach with Internal`Bag

This is very, very slow, also because random numbers a more efficiently created in bulks instead of one-by-one.

n = 1000000;
Do[
   x = RandomReal[{0, 2}];
   y = RandomReal[{0, 1}];
   If[y <= pdf[x], Internal`StuffBag[bag, x]];
   ,
   {n}
   ]; // RepeatedTiming // First
rand = Internal`BagPart[bag, All];

3.2

This takes about 3.2 seconds...

An approach with Compile and Internal`Bag

Compiling the latter can be faster by more than two orders of magnitude, though.

cf = Block[{x},
   With[{pdfx = pdf[x]},
    Compile[{{n, _Integer}},
     Block[{x, y, bag},
      bag = Internal`Bag[Most[{0.}]];
      Do[
       x = RandomReal[{0, 2}];
       y = RandomReal[{0, 1}];
       If[y <= pdfx, Internal`StuffBag[bag, x]];
       ,
       {n}
       ];
      Internal`BagPart[bag, All]
      ],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];

n = 1000000;
rand = Join @@ cf[ConstantArray[n/4, {4}]]; // RepeatedTiming // First

0.022

$\endgroup$
  • $\begingroup$ It's quite faster, but this is the Inversion method, right? I'm trying to learn how to use the ac/rej method $\endgroup$ – Hugo David Costa Pereira Aug 29 '18 at 18:12
  • $\begingroup$ Doh. Okay...... $\endgroup$ – Henrik Schumacher Aug 29 '18 at 18:12
  • $\begingroup$ thanks by the way :) $\endgroup$ – Hugo David Costa Pereira Aug 29 '18 at 18:15
  • $\begingroup$ Ok, i think i got it. But, when it comes to the Internal´Bag section i dont think i understood the true ideia of it. Does it work like the AppendTo function? where you write "bag" it means that this is an empty list where you will attach the tested numbers? $\endgroup$ – Hugo David Costa Pereira Aug 29 '18 at 18:48
  • $\begingroup$ @Hugo Yes, it is exactly as you anticipated. You can search a bit for Internal´Bag on this site. It is not officially documented, though. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 18:49

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.