# Efficient generation of acceptance rejection random number array

Hi there,

I have the following (unnormalized) probability density function (pdf):

f(x) = -x*(x-1)


whose support is the interval (0, 1). I have coded the following to generate a 10000-element array of random numbers drawn from the pdf above:

f_rand = []
for i in 1:10000
x = rand()
u = 1/4*rand()
f_cand = f(x)
if u < f_cand
push!(f_rand, x)
end
end


which does seem to work (I plotted the corresponding histogram(f_rand) for several trials).
I would like to know whether this is coded is an efficient way or I should do some streamling (“vectorization” of the loops, perhaps?).

Here a rejection sampler I once coded: https://github.com/mauro3/KissMCMC.jl/blob/b3c29e1c2c0bfe4ee132bd9640cfe45b33e7f598/src/samplers-serial.jl#L42

Also have a look at the performance tips in the manual. For instance you should put your code into a function, profile it with BenchmarkTools, and check with @code_warntype.

4 Likes

Some things to emulate from @mauro3’s implementation:

• preallocate the result when its size is knowable
• use @inline to modify blocks that index into vectors/arrays
• use broadcasting (the ‘.’ as in a .* vec) rather than
handle items one at a time (unless doing so is intentional)
1 Like

I have tried @mauro3 implementation with:

f_rand = []
for i ∈ 1:10000
push!(z, rejection_sample_unif(x -> x*(1-x), [0.0, 1.0], 1/4.0))
end


and it worked fine.

Following again @mauro3 sugestion, which states it might be faster, I tried using ApproxFun as well, via:

f(x) = -x*(x-1)
f_rand = ApproxFun.sample(f, 10000)


or

f = Fun(x*(1-x))
f_rand = ApproxFun.sample(f, 10000)


Both trials were unsuccessful… I did not manage to understand the use of ApproxFun.sample.

Thanks

I now have made my homework and read the documentation for the package ApproxFun.jl and the working code is:

using LinearAlgebra, Plots, ApproxFun
x = Fun(identity, 0.0..1.0)
f = x*(1 - x)
f_rand = ApproxFun.sample(f, 10000)


The corresponding

histogram(f_rand)


seemed to represent the pdf ok now! Thanks everybody.

1 Like

You can also use the inversion method with F(x) = 6*(x^2/2 - x^3/3)

Fi(y) = real(-1/4*im*(sqrt(3) + -im)*(-2y + 2sqrt(-(y - 1)*y)*im + 1)^(1/3) + (im*(sqrt(3) + im))/(4*(-2y + 2sqrt((1 - y)*y)*im + 1)^(1/3)) + 1/2)

Fi(rand())


Looks like some simplification are possible (it’s going through the complex numbers)

I would start with putting the code in a function and make sure that the array you’re pushing into is typed concretely, right now it’s an array of Any.

4 Likes