# Acceptance Rejection Method

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.

• There is standard answer 2635 to such questions. Commented Jul 10, 2019 at 8:33

# 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 InternalBag

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], InternalStuffBag[bag, x]];
,
{n}
]; // RepeatedTiming // First
rand = InternalBagPart[bag, All];

3.2

## An approach with Compile and InternalBag

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 = InternalBag[Most[{0.}]];
Do[
x = RandomReal[{0, 2}];
y = RandomReal[{0, 1}];
If[y <= pdfx, InternalStuffBag[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

• It's quite faster, but this is the Inversion method, right? I'm trying to learn how to use the ac/rej method
– hcp
Commented Aug 29, 2018 at 18:12
• Doh. Okay...... Commented Aug 29, 2018 at 18:12
• thanks by the way :)
– hcp
Commented Aug 29, 2018 at 18:15
• 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?
– hcp
Commented Aug 29, 2018 at 18:48
• @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. Commented Aug 29, 2018 at 18:49