CODEWITHSUNDEEP
pythonperformancemeanmontecarlo
user2863620
user2863620

Reputation: 643

Calculating mean of sample in monte carlo simulation

I did a Monte Carlo simulation of n samples. For each sample i, I need to calculate the value Xi so probably, the results that I will obtain is:

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

(Here Xi can be a matrix or number).

Now I want to calculate the mean of theses samples that I call Xmean. So I need to obtain something like this:

Xmean = [X1, (X1+X2)/2, (X1+X2+X3)/3 ... , (X1+X2+...+Xn)/n]

In Python, I write a code:

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

It works well but too slow, I would like to know if I can speed up this code? And if you guys could suggest to me some interesting Python's library that help for Monte Carlo simulation.

Thanks,

Upvotes: 0

Views: 2192

Answers (3)

Donbeo
Donbeo

Reputation: 17617

If you use numpy it should be easy.

import numpy as np

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

Upvotes: 0

senshin
senshin

Reputation: 10360

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)))

Result:

mc0                  26.023956370918587
mc1                  0.1423197092108488
mc2                  0.13584513496654083

There's no point in redoing the sum over x(1)..x(i) at each loop iteration, when you already have that information available in xmean. The numpy version by Donbeo is marginally faster than the pure-Python version by me, both of which are nearly 200 times faster (for these data, anyway) than the original version.

Upvotes: 1

Liteye
Liteye

Reputation: 2801

By simply reducing the amount of calculation,

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

Upvotes: 0

Related Questions