Comment on page

# MC Integration

## Classical Monte Carlo Integration

Robert and Casella (2013) introduces the classical Monte Carlo integration:

## Toy Example

$I = \int _D g(\mathbf x)d\mathbf x$

1. 1.
计算区域的volume：
$V = \int_D d\mathbf x$
2. 2.
近似：
$\hat I_m=V\frac{1}{m}\sum\limits_{i=1}^mg(\mathbf x^{(m)})$

$\lim_{m\rightarrow \infty} \hat I_m=I$

$\frac{1}{V}{}\sqrt{m}(\hat I_m-I)\rightarrow N(0, \sigma^2)$

$\sigma^2=var(g(\mathbf x))$

$I = \int_Dg(\mathbf x)d\mathbf x = \int_D Vg(\mathbf x)\frac 1V d\mathbf x = {\mathbb E}_f[Vg(\mathbf x)],$

## 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)