Efficient generation of acceptance rejection random number array

New to Julia
1

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?).

Thanks in advance.

2

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
3

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
  • distinct tasks/task-elements are distinct (smaller) functions
  • use broadcasting (the ‘.’ as in a .* vec) rather than
    handle items one at a time (unless doing so is intentional)
1 Like
4

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

5

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
6

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)

7

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