python - 在python中使用蒙特卡洛方法

Question

我正在处理下图中所写问题的第一个版本。我使用 rand 命令生成了一个随机点，并测试了该点是否在圆圈内。我的代码是否接受蒙特卡洛函数的数量作为输入值？我相信我选择了 N，点数足够小，这样我就不会耗尽内存。此外，当我运行这段代码时，它运行时没有任何错误，但图表没有显示出来。在我可能出错的地方寻求帮助。

import numpy as np
import matplotlib.pyplot as plt
from random import random

xinside = []
yinside = []
xoutside = []
youtside = []

insidecircle = 0
totalpoints = 10**3

for _ in range(totalpoints):
    x = random()
    y = random()
    if x**2+y**2 <= 1:
        insidecircle += 1
        xinside.append(x)
        yinside.append(y)
    else:
        xoutside.append(x)
        youtside.append(y)

fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.scatter(xinside, yinside, color='g', marker='s')
ax.scatter(xoutside, youtside, color='r', marker='s')
fig.show()

Answer 1

没有显示的图表很神秘，也许试试plt.show()。或者，您可以使用 savefig 保存绘图。这是代码第一部分的工作函数(只需修改问题中发布的代码)以及所需的输出图。

import numpy as np
import matplotlib.pyplot as plt
from random import random

def monte_carlo(n_points):
    xin, yin, xout, yout = [[] for _ in range(4)] # Defining all 4 lists together
    insidecircle = 0

    for _ in range(n_points):
        x = random()
        y = random()
        if x**2+y**2 <= 1:
            insidecircle += 1
            xin.append(x)
            yin.append(y)
        else:
            xout.append(x)
            yout.append(y)

    print ("The estimated value of Pi for N = %d is %.4f" %(n_points, 4*insidecircle/n_points))

    fig, ax = plt.subplots()
    ax.set_aspect('equal')
    ax.scatter(xin, yin, color='g', marker='o', s=4)
    ax.scatter(xout, yout, color='r', marker='o', s=4)
    plt.savefig('monte_carlo.png')

n_points = 10**4
monte_carlo(n_points)

> The estimated value of Pi for N = 10000 is 3.1380

矢量化方法 如果您在函数中排除print 语句，您可以将其称为单行代码。我把时间分析留作你的作业

import numpy as np
import matplotlib.pyplot as plt

def monte_carlo(n_points, x, y):
    pi_vec = 4*(x**2 + y**2 <= 1).sum()/n_points
    print ("The estimated value of Pi for N = %d is %.4f" %(n_points, pi_vec))

# Generate points
n_points = 10**4
x = np.random.random(n_points)
y = np.random.random(n_points)
# x = [random() for _ in range(n_points)] # alternative way to define x
# y = [random() for _ in range(n_points)] # alternative way to define y

# Call function
monte_carlo(n_points, x, y)

> The estimated value of Pi for N = 10000 is 3.1284

或者，您可以通过使用单个 x 和 y 点数组来摆脱两个变量 x 和 y，如下所示:

pts = np.random.uniform(0, 1, 2*n_points).reshape((n_points, 2))
monte_carlo(n_points, pts) # Function will be def monte_carlo(n_points, pts):

并在函数中使用

pi_vec = 4*(pts[:,0]**2 + pts[:,1]**2 <= 1).sum()/n_points