初级Python montecarlo模拟

Question

我是Python的初学者，正在练习我们的指导老师安排的练习。我正在努力解决这个问题。在

在Python编辑器中，编写montecarlo模拟来估计数字π的值。 具体来说，请遵循以下步骤： A、 生成两个数组，一个称为x，一个称为y，每个数组包含100个元素， 它们是介于-1和1之间的随机且均匀分布的实数。 B、 以点的形式绘制y和x。相应地标记轴。 C、 写下一个数学表达式来定义哪些（x，y）对数据点 位于半径为1的圆内，以图形的（0，0）原点为中心。 D、 使用布尔掩码来识别圆内的点，并在 在B中已经绘制的数据点上有不同的颜色和标记大小

这就是我现在所拥有的。在

import numpy as np
import math
import matplotlib.pyplot as plt
np.random.seed(12345)
x = np.random.uniform(-1,1,100) 
y = np.random.uniform(-1,1,100) 
plt.plot(x,y) //this works


for i in x:
    newarray = (1>math.sqrt(y[i]*y[i] + x[i]*x[i]))
plt.plot(newarray)

有什么建议吗？在

Answer 1

正如注释中指出的，代码中的错误是for i in x应该是for i in xrange(len(x))

如果你想像语句中所说的那样使用布尔掩码，你可以这样做

    import pandas as pd
    allpoints = pd.DataFrame({'x':x, 'y':y})

    # this is your boolean mask
    mask = pow(allpoints.x, 2) + pow(allpoints.y, 2) < 1
    circlepoints = allpoints[mask]

    plt.scatter(allpoints.x, allpoints.y)
    plt.scatter(circlepoints.x, circlepoints.y)

把点数增加到10000你会得到这样的结果

要估算π，可以使用著名的蒙特卡罗推导法

^{pr2}$

Answer 2

你已经接近解决方案了。稍微改变一下你的形状：

import numpy as np
import math
import matplotlib.pyplot as plt

np.random.seed(12345)

N = 10000
x = np.random.uniform(-1, 1, N) 
y = np.random.uniform(-1, 1, N) 

现在，我们计算一个在这种情况下有意义的标准，例如点到原点的距离：

^{pr2}$

然后我们用Boolean Indexing来区分单位圆内外的点：

q = (d <= 1)

在这一点上，蒙特卡罗假设成立了。我们假设圆和平面上均匀分布点的比率U(-1,1)xU(-1,1)代表单位圆和正方形的面积。然后我们可以从圆/正方形内点的比率来统计地评估pi = 4*(Ac/As)。这导致：

pi = 4*q.sum()/q.size # 3.1464

最后我们绘制结果：

fig, axe = plt.subplots()
axe.plot(x[q], y[q], '.', color='green', label=r'$d \leq 1$')
axe.plot(x[~q], y[~q], '.', color='red', label=r'$d > 1$')
axe.set_aspect('equal')
axe.set_title(r'Monte Carlo: $\pi$ Estimation')
axe.set_xlabel('$x$')
axe.set_ylabel('$y$')
axe.legend(bbox_to_anchor=(1, 1), loc='upper left')
fig.savefig('MonteCarlo.png', dpi=120)

It输出：