Skip to main content
 首页 » 编程设计

python R 实现蒙特卡洛算法计算pi值

2022年07月19日138mengfanrong

算法说明

下面是摘抄内容,说明的非常详细,如果不懂就不要往下看了;

蒙特卡洛算法是通过概率来计算pi的值的。对于一个单位为1的正方形,以其某一个顶点为圆心,边为半径在正方形内画扇形(一个1/4的圆形的扇形),那么扇形的面积就是pi/4。这样,利用概率的方式,“随机”往正方形里面放入一些“”,根据这些点在扇形内的概率(在扇形内的点数/投的总点数),就可以得到扇形的面积。

简单理解,由于正方形为单位1的边长,所以正方形面积S1=1. 假设随机投如N个点,其中M个落入扇形内,那么可以认为扇形的面积为M/N*1。当然,本身,这就是一个概率的过程,不是精确的。当N足够大的时候,得到的面积更精确。由于已经扇形的面积为pi/4,所以根据M/N*1=pi/4就可以计算出pi的值。

这里的关键问题是:(1)如何随机投入点?由于一般语言提供了随机函数,所以利用随机函数就可以模拟随机投点的过程。(2)如何判断点在扇形内?这时候,我们先考虑建立的坐标系,假设单位正方形的顶点在原点,正方形在第一象限,即正方形范围为(0,0)到(1,1)的坐标内,并且假设我们扇形的圆心也为原点,那么根据圆的方程,假设投的点坐标为(x,y),如果x*x+y*y<1,则点(x,y)在扇形内。

参考下图进行理解:

pi = 4*M/N

语言实现

R语言版本

getpi <- function(){ 
 
  N <- 1000000 
  cnt <- 0 
  for (i in 1:N) { 
    x <- runif(1,0,1) 
    y <- runif(1,0,1) 
     
    if((x*x + y*y) <1) 
      cnt= cnt + 1 
     
  } 
   
  return ( 4 * cnt / N) 
} 
 
getpi()


python版本

    def getPi(): 
        cnt = 0 
        for i in range(N) : 
            x = random.uniform(0,1) 
            y = random.uniform(0,1) 
             
            if (x*x + y*y) < 1 : 
                cnt += 1 
        vPi = 4.0 * cnt / N 
        return vPi 
 
    print(getPi())

本文参考链接:https://blog.csdn.net/neweastsun/article/details/51800132