Calculating mean of sample in monte carlo simulation
Question: 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,


Answer: 
    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.

