我正在寻找使用蒙特卡罗模拟的 Pi但不知何故我迷失了自己的代码

问题描述

我正在做一些在线研究和教程，然后当我尝试使用蒙特卡罗找到 pi 的近似值时卡住了。

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

random.seed(1234) 
old_est = 0
n = 1
MC_N= 1000000
results = np.repeat(0,MC_N)

for i in range(MC_N):
    x = random.random() 
    y = random.random()
    A_n = math.sqrt(x*2 + y*2)
    new_est = ((n-1)/n) * old_est + (1/n)*A_n
    results[n] = new_est
    if A_n <= 1:
        n +=1
    if n > MC_N:
        break
    old_est = new_est
    A_n = 4 * n / MC_N
print(results)
plt.plot(results)
plt.ylabel('numbers')
plt.show()

所以我正在使用蒙特卡罗方法，并希望有一个早期停止方法来获得最佳估计。 我想在 3.14 左右，但我得到了 0 有人可以帮助我为什么我错了吗？

PS：它应该是计算pi的近似数，为了估计PI的值，我们需要正方形的面积和圆的面积。为了找到这些区域，我们将在表面上随机放置点并计算落在圆圈内的点和落在正方形内的点。这将为我们提供他们区域的估计数量。因此，我们将使用点数作为面积，而不是使用实际面积。此外，我想表明随着迭代次数的增加，估计误差也呈指数下降。

解决方法

这是一种方法：

import numpy as np
import matplotlib.pyplot as plt
N = 1000000

# N random x,y points drawn uniformly from square(0,1)
x = np.random.rand(N)
y = np.random.rand(N)

# 1 if inside circle radius 1,0 otherwise (first quadrant only)
circle_counts = (x**2 + y**2 < 1).astype(int)

# cumulative fraction of points that fall within circle region
pi_estimate = 4*circle_counts.cumsum()/np.arange(1,N+1)

# pi_estimate = array([4,4,...,3.14215828,3.14215914,3.14216])

plt.plot(x=np.arange(1,N+1),y=pi_estimate)
plt.show()

,

我在维基上查了一下：https://en.wikipedia.org/wiki/Monte_Carlo_method 看来你的算法很不一样，尤其是行new_est = ((n-1)/n) * old_est + (1/n)*A_n。我不确定这是否是您想要的：

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

random.seed(1234) 
MC_N = 100
results = []

est = 0
for i in range(1,MC_N):
    x = random.random() 
    y = random.random()
    A_n = math.sqrt(x**2 + y**2)    # not x*2 + y*2  
    if A_n < 1:
        est += 1    
        
    results.append(4 * est / MC_N)
    
print(results)
plt.plot(results)
plt.ylabel('numbers')
plt.show()