假設有一任意函數 \(g(x)\)，在區間[a,b]中可積分，積分式如下： \[\int^b_ag(x)dx\] 則蒙地卡羅積分可藉由將積分式表達為函數 \(g(x)\)的期望值，再藉由亂數平均近似該期望值，如下式： \[\int^b_ag(x)dx = (b-a)\int^b_ag(x)\frac{1}{b - a}dx = (b-a)E(g(X))\] 其中函數\(g(X)\)的參數 \(X\) 為服從均勻分配的隨機變數。蒙地卡羅積分之所以有效是基於大數法則，當樣本數越大，樣本平均數就會越接近母體分配的期望值。因此根據此法則，生成越多服從均勻分配的亂數代入函數 \(g(X)\)，函數 \(g(X)\) 的平均數將越接近真正的積分值。
若函數為兩個變數的可積分函數，則根據大數法則，可得結果如下： \[ \lim_{n \to \infty}\sum^{n}_{i = 1}g(X_i,Y_i)/n = \int^b_a\int^d_c g(x,y)\frac{1}{(b - a)(d - c)}dxdy \] 其中 \(X_i\)和 \(Y_i\)是獨立的。模擬時先分別生成兩組獨立的亂數， 逐對代入函數 \(g(x)\)後，再求平均數。
蒙地卡羅積分並非只能使用均勻分配，假如 \(X\) 的密度函數為 \(f(x)\)，則函數 \(g(x)\) 的積分可表達為以下形式： \[E(\frac{g(X)}{f(X)}) = \int^b_a\frac{g(x)}{f(x)}f(x)dx =\int^b_ag(x)dx \] 蒙地卡羅積分不是每次都能奏效，如果 \(\frac{g(X)}{f(X)}\) 變化太大可能導致平均數難以收斂。因此決定使用何種 \(f(x)\)將影響是否能求出近似積分值，盡可能選擇能讓 \(\frac{g(X)}{f(X)}\) 保持常數或簡單的函數的分配。
求 \(\int^1_0x^4dx\)
r.unif <- runif(n = 10000)
mean(r.unif^4)*(1 - 0)
## [1] 0.1981422
2 # exact answer
## [1] 2求 \(\int^5_2sin(x)dx\)
r.unif <- runif(n = 10000, min = 2, max = 5)
mean(sin(r.unif))*(5 - 2)
## [1] -0.7109953
-0.7 # exact answer
## [1] -0.7求 \(\int^{10}_3 \int^7_1 sin(x - y) dxdy\)
r.unif.y <- runif(n = 10000000, min = 3, max = 10)
r.unif.x <- runif(n = 10000000, min = 1, max = 7)
mean(sin(r.unif.x - r.unif.y))*(10 - 3)*(7 - 1)
## [1] 0.1014544求 \(\int^{\infty}_{1} exp(-x^2)dx\)
為了調整積分範圍，可令 \(Y = X -1\)，則此積分可表達為\(\int^{\infty}_{0} exp(-(y +1)^2)dy\)
r.exp <- rexp(n = 100000)
mean(exp(-(r.exp + 1)^2) / dexp(x = r.exp))
## [1] 0.1393212\[\int^{\pi}_0 sin(x) dx\]
kN <- 100000
r.unif <- runif(n = kN, min = 0, max = pi)
mean(sin(r.unif))*(pi - 0)
## [1] 2.006404
-(cos(pi) - cos(0)) # true value
## [1] 2\[\int^{\pi}_1 e^x dx\]
kN <- 100000
r.unif <- runif(n = kN, min = 1, max = pi)
mean(exp(r.unif)) * (pi - 1)
## [1] 20.40891
exp(pi) - exp(1) # true value
## [1] 20.42241\[\int^{\infty}_0 e^{-x} dx\]
kN <- 1
r.exp <- rexp(n = kN, rate = 1)
mean(exp(-r.exp) / dexp(x = r.exp, rate = 1))
## [1] 1
-(exp(-Inf) - exp(0)) # true value
## [1] 1\[\int^{\infty}_0 e^{-x^3} dx\]
kN <- 1000000
# 方法1 使用指數分配生成亂數
r.exp <- rexp(n = kN, rate = 1)
mean(exp(-r.exp^3) / dexp(x = r.exp, rate = 1))
## [1] 0.8932502
# 方法2 使用條件標準常態分配生成亂數，給定X>0
r.norm <- rnorm(n = kN)
r.norm <- r.norm[r.norm > 0]
mean(exp(-r.norm^3) / (dnorm(x = r.norm) / 0.5))
## [1] 0.8928141\[\int^{3}_0 sin(e^{x}) dx\]
kN <- 1000000
r.unif <- runif(n = kN, min = 0, max = 3)
mean(sin(exp(r.unif))) * (3 - 0)
## [1] 0.6084642\[\int^{1}_0 \frac{1}{\sqrt{2\pi}} e^{-x^2/2} dx\]
kN <- 1000000
# 方法1 使用均勻分配生成亂數
r.unif <- runif(n = kN, min = 0, max = 1)
r.unif.fun <- exp(-r.unif^2 / 2) * (1 / sqrt(2 * pi))
mean(r.unif.fun)*(1 - 0)
## [1] 0.341374
# 方法2 使用Beta分配生成亂數
r.beta <- rbeta(n = kN, shape1 = 1, shape2 = 1)
r.beta.fun <- exp(-r.beta^2 / 2) * (1 / sqrt(2 * pi))
mean(r.beta.fun / dbeta(x = r.beta, shape1 = 1, shape2 = 1))
## [1] 0.3412591\[\int^{1}_0 \int^{1}_0 cos(x - y) dxdy\]
kN <- 1000000
# 方法1 使用均勻分配生成亂數
r.unif.x <- runif(n = kN)
r.unif.y <- runif(n = kN)
mean(cos(r.unif.x - r.unif.y)) * (1 - 0) * (1 - 0)
## [1] 0.9193045
# 方法2 使用Beta分配生成亂數
r.unif.x <- runif(n = kN)
r.unif.y <- runif(n = kN)
d.fun.x <- dbeta(x = r.unif.x,shape1 = 1, shape2 = 1)
d.fun.y <- dbeta(x = r.unif.y,shape1 = 1, shape2 = 1)
mean(cos(r.unif.x - r.unif.y) / d.fun.x * d.fun.y)
## [1] 0.9194405\[\int^{5}_0 \int^{2}_0 e^{-(x+y)^2} (x+y)^2 dxdy\]
kN <- 1000000
# 方法1 使用均勻分配生成亂數
r.unif.x <- runif(n = kN, min = 0, max = 2)
r.unif.y <- runif(n = kN, min = 0, max = 5)
r.fun <- exp(-(r.unif.x + r.unif.y)^2) * (r.unif.x + r.unif.y)^2
mean(r.fun) * (2 - 0) * (5 - 0)
## [1] 0.4950419
# 方法2 使用條件指數分配生成亂數
# 生成條件指數分配亂數x
r.exp.x <- rexp(n = kN)
r.exp.x <- r.exp.x[r.exp.x < 2]
# 生成條件指數分配亂數y
r.exp.y <- rexp(n = kN)
r.exp.y <- r.exp.y[r.exp.y < 5]
# 將亂數(x,y)代入被積分函數
r.fun <- exp(-(r.exp.x + r.exp.y)^2) * (r.exp.x + r.exp.y)^2
# 分別計算X和Y條件指數分配
d.exp.x.giv <- dexp(x = r.exp.x) / pexp(q = 2)
d.exp.y.giv <- dexp(x = r.exp.y) / pexp(q = 5)
# 求積分值
mean(r.fun / (d.exp.x.giv * d.exp.y.giv))
## [1] 0.494803“A First Course in Statistical Programming with R” [W. John Braun, Duncan J. Murdoch]