R语言蒙特卡罗模拟计算圆周率

发表于

蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。^[蒙特卡罗方法]

1
2
3
4
5
6
7
8
9
10
# x和y构成了一个正方形
# (x, y)表示落在边长为1的正方形的点
x <- runif(10000)
y <- runif(10000)
z <- x^2 + y^2
# sum(z <= 1)统计落在单位圆内的点数
# 由于我们只统计了第一象限的点，但是根据对称性，点落在每个象限的概率是相等的
# 故下面的式子要乘以4
sum(z <= 1) / 10000 * 4
1
#> [1] 3.1472

可视化

1
require(plotrix)
1
#> Loading required package: plotrix
1
require(grid)
1
#> Loading required package: grid
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
simp <- function(n){
  x <- runif(n)
  y <- runif(n)
  xy <- x^2 + y^2
  ind <- xy <=1
  up <- sum(xy <= 1)
  pai <- up / n * 4
  
  # 绘图
  plot(c(0, 1), c(0, 1), type = 'n', asp = 1, 
       main = bquote(pi == .(pai)), xlab="",ylab="")
  rect(0, 0, 1, 1)
  draw.circle(0, 0, 1)
  points(x, y, pch = 20, cex = 0.5,
         col = ind+1)
}
simp(50000)

上面的代码来自豆瓣pi的R语言模拟，有所修改。