## estimate pi
##
## I = \int H(x, y)dxdy
## where H(x, y) = 1 when x^2+y^2 <= 1;
## otherwise H(x, y) = 0
## volume
V = 4
n = 100000
x = runif(n, -1, 1)
y = runif(n, -1, 1)
H = x^2+y^2
H[H<=1] = 1
H[H>1] = 0
I = V* mean(H)
cat("I = ", I, "\n")
## n = 100, I = 2.96
## n = 1000, I = 3.22
## n = 10000, I = 3.1536
## n = 100000, I = 3.14504
Polar Simulation
Note: Unless otherwise stated, the algorithms and the corresponding screenshots are adopted from Robert and Casella (2013).
The following Julia program can be used to do polar simulation.
using Statistics
function polarsim(x::Array)
while true
phi1 = rand()*2*pi
phi2 = rand()*pi - pi/2
u = rand()
xdotxi = x[1]*cos(phi1) + x[2]*sin(phi1)*cos(phi2) + x[3]*sin(phi1)*sin(phi2)
if u <= exp(xdotxi - sum(x.^2)/2)
return(randn() + xdotxi)
end
end
end
# approximate E^pi
function Epi(m, x = [0.1, 1.2, -0.7])
rhos = ones(m)
for i = 1:m
rhos[i] = 1/(2*polarsim(x)^2+3)
end
return(mean(rhos))
end
# example
Epi(10)
