问题描述
我试图从蒙特卡罗方法绘制 pi 的直方图分布,但每次运行模拟时我得到的直方图要么向左倾斜,要么向右倾斜,而不是近似对称的直方图,峰值约为3.14.输出直方图也有一些差距,我认为我正确地逼近了 pi.我的代码如下:
I am trying to plot a histogram distribution of pi from the Monte Carlo method but I am getting histograms that are either skewed to the left or right each time I run the simulation instead of a histogram that is approximately symmetric and peaks at around 3.14. The output histograms also have some gaps and I think I am approximating pi correctly. My code is below:
[...(importing relevant modules)]
N = 1000 #total number of random points
circlex = []
circley = []
squarex = []
squarey = []
pis = []
for i in range(1, M + 1):
x1 = random.uniform(-1,1)
y1 = random.uniform(-1,1)
if x1**2 + y1**2 <=1:
circlex.append(x1)
circley.append(y1)
else:
squarex.append(x1)
squarey.append(y1)
pi = (4.0*len(circlex))/i
pis.append(pi)
print(pi)
print(pis)
plt.hist(pis, color = 'g')
输出:
我遗漏了什么或做错了什么?
What am I missing or doing wrong?
推荐答案
您的代码实际上是正确的.但是有两件事你忘了考虑:
Your code is actually correct. But there are two things you forgot to take into account:
- 绘制直方图可能不是一种有用的可视化方法.实际上,bins 非常大,很难区分正确的近似值(如 3.14)和错误的近似值(如 3.2)
- 这个算法需要很长时间才能收敛到 pi.
作为参考,我使用了相同的方法并得到了那些结果(顺便说一下,您的代码中有一个错字,您应该将 range(1, M + 1)
转换为 范围(1, N + 1)
):
As a reference, I used the same method and get those results (by the way, there's a typo in your code, you should transform range(1, M + 1)
into range(1, N + 1)
):
approx_pi(N=100) # 3.2
approx_pi(N=1_000) # 3.188
approx_pi(N=10_000) # 3.1372
approx_pi(N=100_000) # 3.145
approx_pi(N=1_000_000) # 3.14378
approx_pi(N=10_000_000) # 3.141584
因此,不要害怕为 N
取更大的值以获得更准确的结果.另外,考虑绘制 pi 近似值的演变而不是直方图,以可视化您的近似值是如何收敛的.
Hence, don't be afraid to take larger values for N
to get more accurate results. Plus, consider plotting the evolution of the approximation of pi rather than an histogram to visualize how your approximation is converging.
最后,根据您的目标,可以使用 numpy 获得更快的代码:
Finally, depending on what your goal is, it is possible to get a faster code using numpy:
import numpy as np
pis = 4 * np.cumsum(np.linalg.norm(np.random.random(size=(N, 2)), axis=1)<= 1) / np.arange(1, N + 1)
现在解释:
- 我们定义了一个形状为
(N,2)
的随机数组.这对应于x
和y
的N
个样本.请注意,此处的采样范围在 0 到 1 之间,但这不会改变估计值. - 我们使用
numpy.linalg.norm
和参数axis=1
来计算每个 2 坐标向量的范数. - 如果采样点位于圆内,我们考虑包含
True
的布尔数组,否则考虑False
- 我们对这个数组应用一个累积和.由于在 Python 中,
True
在被视为整数时被视为1
,因此累积和包含在索引i
处的点数仅考虑i
第一个样本时的圆圈. - 我们除以
np.arange(1, N + 1)
,其中包含在索引i
处的相应采样点数. - 我们乘以
4
得到 pi 的近似值.
- We define a random array of shape
(N,2)
. This corresponds to theN
samples ofx
andy
. Note that they are sampled here between 0 and 1, but this does not change the estimation. - We compute the norm of each of these 2-coordinates vectors using
numpy.linalg.norm
with argumentaxis=1
. - We consider the boolean array containing
True
if the sampled point lies within the circle, andFalse
otherwise - We apply a cumulative sum on this array. Since in Python,
True
is considered as1
when considered as an integer, the cumulative sum contains at indexi
the number of points that are in the circle when considering only thei
first samples. - We divide by
np.arange(1, N + 1)
which contains at indexi
the corresponding number of sampled points. - We multiply by
4
to get an approximation of pi.
这段代码真的很丑,但也快得多(比迭代版本快 10 倍左右).我认为根据您的需要,您可能会对此感兴趣.
This code is really ugly, but also way faster (around 10 times faster than the iterative version). I thought that depending on your needs, you may find this of interest.
这篇关于如何绘制蒙特卡罗 pi 直方图?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!