机器学习 numpy求圆周率π（蒙特卡罗法）

解释一下蒙特卡罗法：
蒙特卡罗法的概念可以解释为利用随机试验求解，重复随机取样仿真出复杂的物理和数学系统模型。

欣赏一下这幅图片理解一下蒙特卡罗的思想
以圆为例，，在这个长宽都为1的坐标轴上，四分之一圆的面积S为S =1/4 * π * 1 * 1，π=4 * S，在坐标轴随机出X个点，与原点距离<=1的点（坐标轴中红点）即为圆内的点，个数为t，面积S可视与红点个数除以总个数X。π即等于4 * (t / X)。

用代码实现上述坐标：

import matplotlib.pyplot as plt
import numpy as np
# 运用np.random随机生成1000个点
x=np.random.rand(1000)
y=np.random.rand(1000)
# 转化为坐标格式
X=np.c_[x,y]
# 计算出与原点距离<=1的点
incircle=X[X[:,0]**2+X[:,1]**2<=1]
# 可视化 
plt.figure(figsize=(9,9))
plt.scatter(x,y,c='b',s=8)
plt.scatter(incircle[:,0],incircle[:,1],c='r',s=8)
# 分界线可视化 可注释
t=np.sqrt(1-x**2) 
plt.scatter(x,t,c='k',s=2)
# 求出相似π
PI=4*incircle.shape[0]/X.shape[0]

由于代码只随机出1000个点，计算出的π偏差较大。

—————————————————————————————

x=np.random.rand(100000000)
y=np.random.rand(100000000)

随机10000000个点 得出的π结果近似真实值 (随机数量越大计算量越大，画图耗时时会很长，建议不要超过一千万)

另一种写法：放到一整个圆中
思路：在正方形区域内随机投点，统计落在圆内的点的数目，求出圆内点数与总点数比例求近似π。（画图耗时较长，建议点数不要超过5000）。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# 随机点数
n=1000
# 半径 原点
r=1.0        
a,b=(0.,0.)
# 正方形区域边界
x_min,x_max=a-r,a+r
y_min,y_max=b-r,b+r
# 在正方形区域内随机投点
x=np.random.uniform(x_min,x_max,n)   # 均匀分布
y=np.random.uniform(y_min,y_max,n)   # 均匀分布
# 计算 点到圆心的距离
t=np.sqrt((x-a)**2+(y-b)**2)
# 统计落在圆内的点的数目
index=np.where(t<=r,1,0)
res=sum(index)
# 计算 pi 的近似值
pi=4*res/n
# 可视化
fig=plt.figure(figsize=(9,9))
axes=fig.add_subplot(111) 
for i in range(len(index)):
    if index[i]==1:
        axes.plot(x[i],y[i],'ro',markersize=3,color='r')
    else:
        axes.plot(x[i],y[i],'ro',markersize=3,color='b')
circle=Circle(xy=(a,b),radius=r,alpha=0.1,facecolor='y')
axes.add_patch(circle)
print(pi)