蒙地卡羅法估計圓周率 (Monte Carlo Estimation of pi)

# 蒙地卡羅法
# 圓周率
# 模擬
# Monte Carlo
# pi
# Simulation
# Programming
# plotrix package

分析:

• 本範例使用蒙地卡羅法以估計圓周率 pi。
• 考慮正方形邊長為單位1, 面積為1平方單位. 中間包括圓形, 其半徑為0.5單位，圓形面積 = pi*(0.5^2) = pi/4.
• 在正方形內隨機產生一個點(x, y), 點落在圓型內的機率大約等於扇型面積 / 正方形面積 = pi/4。
• 考慮(x,y)落在圓形中的次數有numberInCircle次, 全部實驗次數為N, 則 pi 大約等於:
  (numberInCircle / N ) * 4.

R程式解說:

• [#3] 首先載入 plotrix 套件, 準備繪製圓形圖.
• [#4] 設定亂數種子, 準備隨機產生資料點(x, y)
• [#5] mfrow=c(1,2) 設定繪圖結果為1列2行
• [#6] 設定 type="n" 表示不繪圖. 因後續須客製化座標軸, axes=FALSE 不繪製座標軸, asp=1設定長寬比例為1. main 標題使用 expression 以標示數學符號pi.
• [#7] 使用 draw.circle {plotrix} 繪製圓形圖.
• [#8] 使用 lines 繪製正外形.
• [#9] 加上4個角落位置的座標軸.
• [#11] 建立 pi.simulation 函數
• [#16-18] 如果點(x, y) 與 (0.5, 0.5) 之距離小於或等於0.5, 表示此點位於圓形之內, 此時 numberInCircle 個數加上1, 繪製紅色點.
• [#21] 點(x, y)不位於圓形內, 繪製藍色點.
• [#23] 計算估計pi值並儲存於 c 物件.
• [#25] 回傳估計pi值物件c.
• [#30] 取出最後一次計算之估計pi值.
• [#39] 使用 par 函數將繪圖結果還原成1列,1行.

蒙地卡羅法估計圓周率圖:

R程式:

# title: Monte Carlo Estiamtion of pi
# date: 2017.10.9
library(plotrix)
set.seed(123)
op <- par(mfrow=c(1,2))
plot(0.5, 0.5, xlim=c(0, 1), ylim=c(0,1), type="n", axes=FALSE, asp=1, xlab="", ylab="", main=expression(paste("Monte Carlo for ", pi)))
draw.circle(0.5, 0.5, radius=0.5)
lines(x=c(0,1,1,0,0), y=c(0,0,1,1,0))
text(c(0.05,0.95,0.95,0.05), c(0.05,0.05,0.95,0.95), c("(0,0)", "(1,0)", "(1,1)", "(0,1)"))

pi.simulation <- function(samplesize) {
 c <- rep(0, samplesize)
 numberInCircle <- 0
 for (i in 1:samplesize) {
 x <- runif(2)
 if (sqrt((x[1]-0.5)^2 + (x[2]-0.5)^2) <= 0.5) {
 numberInCircle <- numberInCircle + 1
 points(x[1], x[2], col="red", pch=".")
 }
 else {
 points(x[1], x[2], col="blue", pch=".")
 }
 c[i] <- (numberInCircle / i) * 4
 }
 return(c)
}

size <- 10000
pi.sim <- pi.simulation(size)
estimation.pi <- pi.sim[size]
text(0.5, 0.5, paste0("sample size=", size), cex=1.5)
text(0.5, 0.4, paste("pi=",estimation.pi))
plot(pi.sim, type="l", 
 main="Monte Carlo method for pi",
 xlab="samples", ylab=expression(paste("Estimation of ", pi)), 
 ylim=c(0,6))
abline(h=estimation.pi, col="red", lty=3)
grid()
par(op)
# end