蒙特卡罗方法在Python中的应用

Question

我一直在尝试使用Python创建一个脚本，让我可以生成大量的点来用于Monte Carlo方法计算估计值Pi。目前我拥有的脚本如下：

import math
import random
random.seed()

n = 10000

for i in range(n):
    x = random.random()
    y = random.random()
    z = (x,y)

    if x**2+y**2 <= 1:
        print z
    else:
        del z

到目前为止，我已经能够生成所有所需的点，但我想要得到的是在运行脚本时产生的点数，以便进行后续计算。我不需要非常精确的结果，只需要一个足够好的估计值。如果有任何建议，将不胜感激。

- Matt Johnson

---

你想要计算圆内有多少个随机对吗？如果是这样，只需使用一个计数器... - zenpoy

Answer 1

如果你正在进行任何形式的重型数字计算，请考虑学习numpy。使用numpy设置，你的问题基本上是一个一维线性问题：

import numpy as np

N   = 10000
pts = np.random.random((N,2))

# Select the points according to your condition
idx = (pts**2).sum(axis=1)  < 1.0
print pts[idx], idx.sum()

给予：

[[ 0.61255615  0.44319463]
 [ 0.48214768  0.69960483]
 [ 0.04735956  0.18509277]
 ..., 
 [ 0.37543094  0.2858077 ]
 [ 0.43304577  0.45903071]
 [ 0.30838206  0.45977162]], 7854

最后一个数字是计算的事件数量，即半径小于一的点的数量。

- Hooked

Answer 2

不确定这是否符合您的要求，但是您可以在range上运行enumerate并获取迭代中的位置：

In [1]: for index, i in enumerate(xrange(10, 15)):
   ...:     print index + 1, i
   ...:
   ...:
1 10
2 11
3 12
4 13
5 14

在这种情况下，index + 1代表当前正在创建的点（index本身是给定迭代开始时创建的点的总数）。此外，如果您使用的是Python 2.x，则xrange通常更适合这些迭代，因为它不会将整个列表加载到内存中，而是根据需要访问它。

- RocketDonkey

---

@ahans 没问题 - 我记得第一次看到它时就想着“啊，原来有这样的方法。”使用它时玩得开心！ - RocketDonkey

Answer 3

只需在循环之前添加hits变量，将其初始化为0，并在if语句内递增hits。最后，您可以使用hits和n计算PI值。

import math
import random
random.seed()

n = 10000
hits = 0  # initialize hits with 0

for i in range(n):
    x = random.random()
    y = random.random()
    z = (x,y)

    if x**2+y**2 <= 1:
        hits += 1
    else:
        del z

# use hits and n to compute PI

- pm007