Acceptance Rejection Method

Question

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.

Answer 1

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]

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 they will have quadratic complexity and be 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