using Plots

Monte Carlo Method

sampling an arbitrary distributions $\rho(x)$

Von Neumann Rejection

generate a pair of random numbers (x, y)

if y < rho(x), accepts x

Direct Sampling

construct $r(x) = F(x) = \int^x dx^\prime \, \rho(x^\prime)$ such that

$ dr = dx \, \rho(x) $;

invert to get $x(r)$

# rejection

function reject(rho; Nx=8000, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0)
    
    i1 = 0
    bag = []
    while true
        
        xx = xmin + rand()*(xmax-xmin)
        yy = ymin + rand()*(ymax-ymin)
        
        if yy < rho(xx)
            push!(bag, xx)
            i1 += 1
        end
        if i1 > Nx
            break
        end
    end
    return bag
end

reject (generic function with 1 method)

# random distribution by direct sampling

tau = 1.0
r2x = r -> -tau *log(1-r)

N = 25000
r_arr = rand(N)

# a new random distribution is just as easy as this:
x_arr = r2x.(r_arr)

x1_arr = reject(x->exp(-x/tau)/tau, xmax=4.0)

histogram(x_arr,normalize=:pdf, label="direct sampling")
plot!(range(0.0, 10.0, length=50), t->exp(-t/tau)/tau, xlim=(0,5), label="analytic")
histogram!(x1_arr, normalize=:pdf, label="rejection")

# application

tau = 1.0
r2x = r -> -tau *log(1-r)

N = 2000
r_arr = rand(N)

# as easy as that
x_arr = r2x.(r_arr)

rho = x -> exp(-x/tau)/tau
f = x -> exp(-x^2)
ret1 = sum( @. f(x_arr)/rho(x_arr) )/N

println([ret1, sqrt(pi)/2])

[0.8761322034645804, 0.8862269254527579]

# random distribution by direct sampling

r2x = r -> tan(pi *(r-0.4))

N = 1500
r_arr = rand(N)

# as easy as that
x_arr = r2x.(r_arr)

histogram(x_arr, normalize=:pdf, label="direct sampling")
plot!(range(-10.0, 10.0, length=2500), x->1/(1+x^2)/pi, xlim=(-10,10), label="analytic")

# application

r2x = r -> tan(pi *(r-0.5))

N = 25000
r_arr = rand(N)

# as easy as that
x_arr = r2x.(r_arr)

rho = x -> 1/(1+x^2)/pi
f = x -> exp(-x^2)
ret1 = sum( @. f(x_arr)/rho(x_arr) )/N

println([ret1, sqrt(pi)])

[1.7799889262259507, 1.7724538509055159]