R语言入门之 – 如何通过Monte Carlo来计算 PI?

JustYY.com 小赖子的英国生活和资讯

8 years ago

上次开始步入R语言的世界, 感觉R还是挺简洁强大的. 学一门程序最好的办法就是敲代码, 敲例子. 在工作生活中如果遇到需要敲代码的时候就得问问自己能否拿R语言来解决? 这样能更好的进步.

我们都知道圆周率可以通过随机在一个正方形(坐标X/Y均为0到1)撒足够多的点. 统计一下点在1/4的圆内(半径为1)的个数和总的撒点个数这个值就会很接近 . 因为圆的面积公式为

这种方法也称之为蒙特卡罗方法, 是一种随机, 统计的方法.

我们可以通过 runif 来生成随机的点, 参数指定点的个数,

x=runif(100000)
y=runif(100000)

每个随机值是在0到1之间的浮点数(也可以指定 min=0,max=1). 然后可以把长度放在另一个向量里:

z=sqrt(x^2+y^2)

这时候我们只要统计出这个z数组里小于或等于1的个数即可. R语言里的 which 函数返回了数组里满足条件的索引值, 那么我们只需要通过:

length(which(z<=1))

即可以得到在圆内点的个数, 记得乘于4再除于样本的个数就可以得到圆周率的近似值 (点数越多精度越接近, 但是个无理数)

length(which(z<=1))*4/length(z)
[1] 3.1454

把在圆内的点画出来 (不指定颜色默认为黑色)

plot(x[which(z<=1)],y[which(z<=1)],xlab="X",ylab="Y",main="Monte Carlo")

再把剩下的点用蓝色画出来 (注意这里用的是 points 方法在原来的图上画):

points(x[which(z>1)],y[which(z>1)],col='blue')

最后面就很直观的能解释这种算法.

monte-carlo-pi-in-r

这个例子中R非常强大, 没用到 for/while (Python, Matlab 也可以). 传统编程语言基本都是需要用到for或者while的. 从这个例子中是不是能看出统计学出生的R语言的一些编程风格呢?

微信公众号: 小赖子的英国生活和资讯 JustYYUK