# Efficient generation of acceptance rejection random number array
**URL:** https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374
**Category:** New to Julia
**Created:** [November 18, 2020, 12:07pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374 "2020-11-18T12:07:49Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@mocalvao](https://discourse.julialang.org/u/mocalvao)
#### Post date: [November 18, 2020, 12:07pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/1 "2020-11-18T12:07:49Z")
Hi there,
I have the following (unnormalized) probability density function (pdf):
```julia
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:
```julia
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.
---
### Author:  [@mauro3](https://discourse.julialang.org/u/mauro3)
#### Post date: [November 18, 2020, 12:24pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/2 "2020-11-18T12:24:42Z")
Here a rejection sampler I once coded: [https://github.com/mauro3/KissMCMC.jl/blob/b3c29e1c2c0bfe4ee132bd9640cfe45b33e7f598/src/samplers-serial.jl#L42](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`.
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [November 18, 2020, 12:47pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/3 "2020-11-18T12:47:43Z")
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)
---
### Author:  [@mocalvao](https://discourse.julialang.org/u/mocalvao)
#### Post date: [November 18, 2020, 1:34pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/4 "2020-11-18T13:34:31Z")
I have tried @mauro3 implementation with:
```julia
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:
```julia
f(x) = -x*(x-1)
f_rand = ApproxFun.sample(f, 10000)
```
or
```julia
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
---
### Author:  [@mocalvao](https://discourse.julialang.org/u/mocalvao)
#### Post date: [November 18, 2020, 2:42pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/5 "2020-11-18T14:42:17Z")
I now have made my homework and read the documentation for the package ApproxFun.jl and the working code is:
```julia
using LinearAlgebra, Plots, ApproxFun
x = Fun(identity, 0.0..1.0)
f = x*(1 - x)
f_rand = ApproxFun.sample(f, 10000)
```
The corresponding
```julia
histogram(f_rand)
```
seemed to represent the pdf ok now! Thanks everybody.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [November 18, 2020, 5:51pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/6 "2020-11-18T17:51:19Z")
You can also use the inversion method with `F(x) = 6*(x^2/2 - x^3/3)`
```julia
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)
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [November 18, 2020, 6:06pm UTC](https://discourse.julialang.org/t/efficient-generation-of-acceptance-rejection-random-number-array/50374/7 "2020-11-18T18:06:04Z")
> [@mocalvao](#):
>
> I would like to know whether this is coded is an efficient way or I should do some streamling (“vectorization” of the loops, perhaps?).
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.