初学者 Python Monte Carlo 模拟

Question

我是 Python 的初学者，正在完成我们的讲师设置的练习。我正在为这个问题而苦苦挣扎。

在Python编辑器中，编写一个Monte Carlo模拟来估计数字π的值。具体执行以下步骤： 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 你会得到这样的东西

要估计圆周率，您可以使用著名的蒙特卡洛推导

    >>> n = 10000
    >>> ( len(circlepoints) * 4 ) / float(n)
    <<< 3.1464

Answer 2

您已接近解决方案。我稍微改造了你的 MCVE:

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)

现在，我们计算一个在此上下文中有意义的标准，例如点到原点的距离：

d = x**2 + y**2

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

q = (d <= 1)

此时 Monte Carlo 假设成立。我们假设圆和平面中均匀分布点的比例 U(-1,1)xU(-1,1) 代表单位圆和正方形的面积。然后我们可以从Circle/Square内的点数比例统计评估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)

它输出：

初学者 Python Monte Carlo 模拟

Beginner Python Monte Carlo Simulation

python

montecarlo