早上中级微观经济学课上复习泰勒展开和麦克劳林展开，顺带讲到了用蒙特卡洛方法实现计算π值，于是下午着手用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()
