Python实现利用蒙特卡洛法求Pi的值

Posted on | Edited on | In | Comments:
| Symbols count in article: 8k | Reading time ≈ 7 mins.

蒙特卡洛法的介绍：

蒙特卡洛（Monte Carlo）法是一类随机算法的统称。

蒙特卡罗方法解题过程的三个主要步骤：

（1）构造或描述概率过程

对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过 程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。

（2）实现从已知概率分布抽样

构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。

（3）建立各种估计量

一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。

计算pi的n种方式

在近代，pi的计算可以巧妙的利用反三角函数和级数积分进行计算，这便是莱布尼茨级数的收敛结果，可以得到

{% raw %}

当然，更好玩的方式是通过频率-概率计算，那么如何实现呢？

很显然在连续的几何概型中：

如果一个单位圆在坐标轴的正中央，被一个2x2的正方形刚好包裹，那么在上面任意投掷一点，在圆内的概率便是pi/4,很显然，这类概率可以通过频率去逼近。

对于大量，重复的随机点分布，可以利用Python的Turtle库和Random库，都是非常基础的Python包，那么代码实现如下：

#Radius is half of the sqare's length
import turtle
import math
import random as r

a = turtle.Turtle()

a.color('blue')

def pink():
    color = (1,r.random(),1)
    return color

def circle():
    '''circle drawing'''
    a.penup()
    a.goto(50,0) #circle center align
    a.pendown()
    a.circle(50) #circle radius
    a.begin_fill()
    a.fillcolor(pink())
    a.end_fill()

b = turtle.Turtle()

b.color('red')

def square():
    '''square drawing'''
    for i in range(4):
        b.forward(100) #square length
        b.left(90)
        b.begin_fill()
        b.fillcolor(pink())
        b.end_fill()

c = turtle.Turtle()

def blue():
    c.begin_fill()
    c.fillcolor('blue')
    c.end_fill()
    c.penup()


def red():
    c.begin_fill()
    c.fillcolor('red')
    c.end_fill()
    c.penup()


def Montecarlo(n):
    '''Points Generate'''
    points = []
    count = 0
    for i in range(n):
        raw = r.random()*100
        colomn = r.random()*100
        xy = (raw,colomn)
        points.append(xy)
        #print(xy)
        for x,y in points:
            distance = math.sqrt((x-50)**2+(y-50)**2)
        if distance <= 50: #in circle
            blue()
            c.goto(x,y)
            c.stamp()
            c.goto(x,y)
            c.up()
            count += 1
        else: #out of circle
            red()
            c.goto(x,y)
            c.stamp()
            c.goto(x,y)
            c.up()
    PI = count*4/n
    return PI

if __name__ == '__main__':
    wn = turtle.Screen()
    wn.tracer(12,25)
    wn.setworldcoordinates(0,0,100,100)
    circle()
    square()
    print("Pi approximate to:" + str(Montecarlo(2000)))

终端输出结果为

Pi approximate to:3.1468
-------------EndThanks for reading-------------
0%