MC Integration

Classical Monte Carlo Integration

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

Toy Example

为了计算积分
我们通常采用MC模拟:
  1. 1.
    计算区域的volume:
  2. 2.
    近似:
根据大数律有
并且由中心极限定理有
其中
另外,注意到
I=Dg(x)dx=DVg(x)1Vdx=Ef[Vg(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)],
则可以直接利用上一节的结论。
1
## estimate pi
2
##
3
## I = \int H(x, y)dxdy
4
## where H(x, y) = 1 when x^2+y^2 <= 1;
5
## otherwise H(x, y) = 0
6
7
## volume
8
V = 4
9
10
n = 100000
11
x = runif(n, -1, 1)
12
y = runif(n, -1, 1)
13
H = x^2+y^2
14
H[H<=1] = 1
15
H[H>1] = 0
16
I = V* mean(H)
17
cat("I = ", I, "\n")
18
19
## n = 100, I = 2.96
20
## n = 1000, I = 3.22
21
## n = 10000, I = 3.1536
22
## n = 100000, I = 3.14504
Copied!

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.
1
using Statistics
2
3
function polarsim(x::Array)
4
while true
5
phi1 = rand()*2*pi
6
phi2 = rand()*pi - pi/2
7
u = rand()
8
xdotxi = x[1]*cos(phi1) + x[2]*sin(phi1)*cos(phi2) + x[3]*sin(phi1)*sin(phi2)
9
if u <= exp(xdotxi - sum(x.^2)/2)
10
return(randn() + xdotxi)
11
end
12
end
13
end
14
15
# approximate E^pi
16
function Epi(m, x = [0.1, 1.2, -0.7])
17
rhos = ones(m)
18
for i = 1:m
19
rhos[i] = 1/(2*polarsim(x)^2+3)
20
end
21
return(mean(rhos))
22
end
23
24
# example
25
Epi(10)
Copied!