MC Integration
Classical Monte Carlo Integration
Toy Example
为了计算积分 d\mathbf%20x)
我们通常采用MC模拟:
- 1.计算区域的volume：
- 2.近似：
根据大数律有
并且由中心极限定理有
其中))
另外，注意到
则可以直接利用上一节的结论。
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!
Last modified 3yr ago
Copy link