用蒙特卡洛方法算pi-基于python和R语言
最近follow了MOOC上一门python课,开始学Python。同时,买来了概率论与数理统计,准备自学一下统计。(因为被鄙视过不是统计专业却想搞数据分析)
有趣的是书里面有一块讲蒲丰投针计算Pi,这是一种随机模拟法,也就是蒙特卡洛法。蒲丰投针之于我太难,暂时没想到怎么用计算机模拟这一过程。
python课中,老师也提到用随机模拟法,也就是蒙特卡洛法(MonteCarlo),用计算机模拟几千次实验,计算pi的近似值。好巧。
就拿python课中的方法,来近似计算pi,分别用python和R实现一下。
至于实验是怎样的,截图老师的PPT。
我用的是python 3
python代码:
from random import random from math import sqrt from time import clock darts=2**22 hist=0 clock() for i in range(1,darts): x,y=random(),random() dist=sqrt(x**2+y**2) if dist<=1.0: hist=hist+1 pi=4*(hist/darts) print('pi is %s'%pi) print('elaspe is %ss'%clock())
python运行结果:
pi is 3.143444061279297 elaspe is 85.991785
R 代码:
#蒙特卡洛方法求pi hist <- 0 darts <- 2^22 start <-proc.time() for (i in 1:darts){ x <- runif(1);y <- runif(1) if(sqrt(x^2+y^2)>1) next hist=hist+1 } pi <- 4*hist/darts proc.time()-start print(paste0('pi is ',pi))
R运行结果:
> proc.time()-start 用户 系统 流逝 31.537 2.477 34.153 > print(paste0('pi is ',pi)) [1] "pi is 3.14076137542725"
总结:R和python都挺好用的,下一步准备试着用python写点小爬虫程序。