怎样用python计算π的值_IV.python初探日记：python实现蒙特卡洛方法计算π值

早上中级微观经济学课上复习泰勒展开和麦克劳林展开，顺带讲到了用蒙特卡洛方法实现计算π值，于是下午着手用python尝试着实现了一下，并用matplotlib输出了一部分数据。

完整的代码在文末，本文适合小白看，完全白纸的都可以，也希望大神们不吝赐教。

一、最简单的实现方法

下面是最简单的实现方式，模拟试验一千万次，但模拟出来的π值并不精确。

import random

zongshu = 10000000

jishu = 0

for i in range(zongshu) :

x = random.random()

y = random.random()

if (x ** 2 + y ** 2) < 1 :

jishu+=1

print(jishu * 4.0 / zongshu)

二、关于代码的效率问题

上面这段代码的计算速度和精度都相当粗糙，然而我还没有对改进代码的想法，所以下面就直接暴力模拟了。下面是粗略地试验了一下各种实验次数下的实现时间，此处引用的是我修改后的代码，原代码为——使用Python语言的蒙特卡洛方法计算圆周率π的一种实现。

from __future__ import division

import random

import time

num = 1

for j in range(1, 7):

startT = time.clock()

# 落入圆内计数

counter = 0

# 往正方形中扔了 10 的 j 次方次

for i in range(10 ** j):

x = random.uniform(-1, 1)

y = random.uniform(-1, 1)

# 落入了圆内

if x**2 + y**2 < 1:

counter = counter + 1

endT = time.clock()

print ( num)

print ( 'pi:{0}'.format(4 * (counter / 10 ** j)))

print ( 'Time:{0}'.format(endT - startT))

print ('')

num += 1

模拟试验一千万次以上的时间就相当漫长了。

三、matplotlib输出左图输出的是模拟的1/4圆试验，右图输出的是试验次数和计算π值的关系，这里主要碰到了下面这几个问题。

1、生成随机点

主要的思想是生成两组随机数分别加入两组数组。原本不知道matplotlib是可以直接输出数组的，折腾了好久才发现这么简单......尝试过的方法有：循环“生成一个随机点”这个步骤；生成两对随机数组，然后一一匹配得到随机点；还有一些奇葩的就不列举了......总之这些步骤都没有去搜索，纯粹想锻炼一下自己的思维。

2、1/4圆“那条线”

一开始是想用matplotlib直接画一个圆，结果发现matplotlib不能直接画圆(也有可能是我没找到，知道的同学告诉我一声)；随后想到圆环、实心圆涂色等几个方法，但实现起来都较为复杂，最后我找到了这样一个方法numpy.meshgrid，直接把这个1/4圆描了出来：

x = y = np.arange(-1, 1, 0.001)

x, y = np.meshgrid(x,y)

p1.contour(x, y, x**2 + y**2, [1])

3、右图代码的取舍

右边的关系图，考虑到速度和精度，取的是1000。

C = [ (j+1) for j in range(1000)]

D = [ M((j+1)) for j in range(1000)]

其实可以改进为下面这样，但我原代码中没有修改。

C = [ 100 * (j+1) for j in range(100)]

D = [ M(100 * (j+1)) for j in range(100)]

4、排版问题

用的是“121”的排版方式，出来的图手工拉动了一下大小，不过似乎matplotlib是可以直接控制输出的，这点我要再琢磨下。

5、遗留问题

左图其实可以直接生成随机点；“1/4圆”的问题没有从根本上解决；输出的排版没有解决；代码的优化完全没有涉及；另外比较现实的是，如何把“1/4圆”一下的点抹成红色，“1/4圆”以上的点保持蓝色，这点我还没有去尝试解决。上面的遗留问题如果有进一步的探索的话，我会在这里跟进的。

四、完整代码

最后是完整的代码。

import random

import matplotlib.pyplot as plt

import numpy as np

A = []

B = []

cishu = 100

for i in range(cishu):

X = random.random()

A.append(X)

Y = random.random()

B.append(Y)

def M(zongshu):

jishu = 0.0

for j in range(zongshu):

x = random.random()

y = random.random()

if ( x**2 + y**2 ) < 1:

jishu += 1

return (jishu/zongshu * 4.0)

fig=plt.figure()

p1=fig.add_subplot(121)

p1.axis([0,1,0,1])

p1.scatter(A, B)

x = y = np.arange(-1, 1, 0.001)

x, y = np.meshgrid(x,y)

p1.contour(x, y, x**2 + y**2, [1])

p2=fig.add_subplot(122)

C = [ (j+1) for j in range(1000)]

D = [ M((j+1)) for j in range(1000)]

p2.plot(C, D)

plt.show()