在蒙特卡罗模拟中计算样本的平均值

Question

我对n个样本进行了蒙特卡罗模拟。对于每个样本 i，我需要计算值 Xi，因此，我将获得的结果很可能是：

X = [X1， X2， ...， Xn]

（这里的Xi可以是矩阵或数字）。

现在我想计算这些样本的平均值，我称之为 Xmean。所以我需要获得这样的东西：

Xmean = [X1， （X1+X2）/2， （X1+X2+

X3）/3 ... ， （X1+X2+...+Xn）/n]

在Python中，我编写了一段代码：

for i in range(N):
    for j in range(i+1):
         Xmean(i) = Xmean(i) + X(j)
    Xmean(i) = Xmean(i) / (i+1)

它运行良好但太慢，我想知道我是否可以加快这段代码？如果你们能向我推荐一些有趣的 Python 库，帮助蒙特卡罗模拟。

谢谢

Answer 1

import timeit, numpy
setup = '''
from __main__ import mc0, mc1, mc2
import random, numpy
random.seed(0)
n = 10**3
data = [random.randint(0, 2**32-1) for _ in range(n)]
np_data = numpy.array([float(x) for x in data])
'''
# your implementation
def mc0(data):
    xmean = []
    for i in range(len(data)):
        xmean.append(0)
        for j in range(i+1):
            xmean[i] += data[j]
        xmean[i] = xmean[i] / (i+1)
    return xmean
# my implementation
def mc1(data):
    xmean = []
    for i, x in enumerate(data):
        if i == 0:
            new = x
        else:
            new = x/(i+1) + xmean[i-1] * (i/(i+1))
        xmean.append(new)
    return xmean
# Donbeo's numpy implementation
def mc2(data):
    xmean = numpy.cumsum(data) / numpy.array(range(1, len(data)+1))
    return xmean

number = 100
things = [('mc0', 'mc0(data)'),
          ('mc1', 'mc1(data)'),
          ('mc2', 'mc2(np_data)')]
for note, call in things:
    print('{:20} {}'.format(note,
                            timeit.timeit(call, setup=setup, number=number)))

结果：

mc0                  26.023956370918587
mc1                  0.1423197092108488
mc2                  0.13584513496654083

在

每次循环迭代中重做x(1)..x(i)的总和是没有意义的，因为您已经在 xmean 中提供了该信息。Donbeo 的 numpy 版本比我的纯 Python 版本略快，两者都比原始版本快近 200 倍（无论如何，对于这些数据）。

Answer 2

通过简单地减少计算量，

Xmean(0) = X(0)
for i in range(N-1):
    Xmean(i+1) = (Xmean(i)*(i+1) + X(i+1))/(i+2)

Answer 3

如果你使用numpy，它应该很容易。

import numpy as np
X = [1,5,3,8,6,9]
Xmean = np.cumsum(X)
Xmean = Xmean/np.array(range(1,len(X)+1)