計算蒙特卡洛模擬中樣本的平均值

我做了n個樣本的蒙特卡洛模擬。對於每個樣品I，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庫來幫助蒙特卡羅模擬。

感謝，

來源

2014-01-18 user2863620

通過簡單地減少了計算量，

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

來源

2014-01-18 12:31:44 Liteye

謝謝，它確實更快：D – user2863620

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的版本比我剛剛發佈的純Python版本要快一些（對於這些數據，無論如何）比原始版本要快200倍。

來源

2014-01-18 12:34:06 senshin

感謝您的回答！ – user2863620

如果你使用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)

來源

2014-01-18 12:56:22 Donbeo

這會引起一個零除錯誤。 – senshin

我認爲這是正確的： Xmean = Xmean/np.array（範圍（1，len（X）+1）） – user2863620

你是對的我改變了代碼 – Donbeo

計算蒙特卡洛模擬中樣本的平均值

回答

相關問題