How to simulate this Gamma expansion in a Python way

Question

Here is the simulation that I want to do:

1. For each of the 10 million simulation paths,
2. I have n = 100 lambda values in sequence (the lambda vector is the same for all paths),
3. Using each of the lambda values, I need to generate a random number from Poisson distribution with mean equals lambda;
4. Say the random numbers from Poisson distribution are [Poi(1), Poi(2), Poi(3)... Poi(100)],
5. Then I need to generate 100 set of random numbers from exponential distribution. In each set, Poi(i) is the number of random number I need.
6. Say the set of random numbers are: [Exp(1,1), Exp(1,2), ..., Exp(1, Poi(1)) Exp(2,1), Exp(2,2), ..., Exp(2, Poi(2)) ... Exp(100,1), ..., Exp(100,Poi(100))]
7. The result that I need for each simulation is one number, which is a weighted average of the sum of each set of exponential random numbers. i.e. result = w(1)*[Exp(1,1)+Exp(1,2)+...+Exp(1,Poi(1))] +...+ w(100)[Exp(100,1)+Exp(100,2)+...+Exp(100,Poi(100))]

How can I do it in a Python way? I write the following code which manage to do one simulation in each for loop, it takes 1 sec for each simulation. Is there a faster way for the simulation?

import numpy as np 
from math import pi 
import time

start = time.clock()

kau = 6.21 
theta = 0.019 
sigma_v = 0.61 
t = 1 
nSim = 10000000 
result = []

n = np.arange(1,7) 
lamb = 16 * pi**2 * n**2 / sigma_v**2 / t / (kau**2 * t**2 + 4 * pi**2 * n**2) 
gamma = (kau**2 * t**2 + 4 * pi**2 * n**2) / 2 / sigma_v **2 / t**2

np.random.seed(1)

for sim in range(nSim):
    Nn = np.random.poisson(lam=lamb)
    y = [np.random.exponential(1,size=i) for i in Nn]
    z = [sum(y[i]) for i in range(len(y))]
    X1 = sum(z/gamma)
    result.append(X1)

print("time elapsed = ",time.clock()-start)

Answer 1

Ill have a crack at this later if I have time but basically the best thing to do in python is use all of your RAM for example if you want to generate 10mm RVs DO NOT do this:

for i in range(10e6):
    rv = stats.norm.rvs(size=1)

instead DO this:

rvs = stats.norm.rvs(size=10e6)

avoiding python loops is much more preferable, so having the variables in existence before use is better (if RAM size permits). If your equations/mechanic permits it you can vectorize all calculations so that they all operate at once and return an array of your desired values.

You have these loops:

for sim in range(nSim) and y = [np.random.exponential(1,size=i) for i in Nn] z = [sum(y[i]) for i in range(len(y))], which are likely to be much faster if you use index slicing of appropriate variables than the inherent loops (e.g. even if you over generate random variables and store in memory and subselect from those it will be quicker than generating precisely those you need in loops).

Edit: for steps 2-7

Here is an example of how to create a mask to subselect the top $n_i$ elements of a random array and sum them:

>>> np.random.seed(1)
>>> rvs = np.random.rand(16).reshape(4,4)
array([[4.17022005e-01, 7.20324493e-01, 1.14374817e-04, 3.02332573e-01],
       [1.46755891e-01, 9.23385948e-02, 1.86260211e-01, 3.45560727e-01],
       [3.96767474e-01, 5.38816734e-01, 4.19194514e-01, 6.85219500e-01],
       [2.04452250e-01, 8.78117436e-01, 2.73875932e-02, 6.70467510e-01]])
>>> index = np.tile(np.arange(1, 5)[:, np.newaxis], (1,4))
array([[1., 1., 1., 1.],
       [2., 2., 2., 2.],
       [3., 3., 3., 3.],
       [4., 4., 4., 4.]])
>>> number_in_col = np.array([1,2,1,3])
>>> mask = number_in_col >= index
array([[ True,  True,  True,  True],
       [False,  True, False,  True],
       [False, False, False,  True],
       [False, False, False, False]])
>>> e_sum = np.sum(rvs * mask, axis=0)
array([4.17022005e-01, 8.12663088e-01, 1.14374817e-04, 1.33311280e+00])
>>> weights = np.array([1.1, 1.2, 1.3, 2.2])
>>> simulation_value = np.sum(e_sum * weights)
4.36691675844615

The above code is essentially the skeleton for steps 2-7.

To do step 1, i.e repeat this 10mm times, you should consider the RAM usage and add in a third axis to your arrays where the third axis contains the information about each simulation. Instead of returning a single simulation value 8.112136084 you will get a 1D array of simulation values dependent upon how many you can muster in one pass.

For example, having memory for index(8-byte) rvs(8-byte) and mask(1-byte) with 100 columns and I'm guessing 50 rows accounts for around (17*100*50) 85,000 bytes per simulation. If you have at most 10GB RAM to spare you could do 100,000 simulations in one go and write the data, then loop through 100 times:

output = np.empty(shape=(100,100000))    
for i in range(100):
    # edit the above for third axis with 100,000 simulations each run
    output[i, :] = simulation_value[:]
return output.reshape(-1,)  # <- output is a 1D array of 10mm values.

I'm not going to test it but I guarantee it will be a lot faster than your current implementation.

Edit to extend 3rd axis for help with step 1

>>> rvs_ext = np.tile(rvs[:,:,np.newaxis], (1,1,3)) # <- create 3 copies along 3rd axis
>>> mask_ext = np.tile(mask[:,:,np.newaxis], (1,1,3))
>>> e_sum_ext = np.sum(rvs_ext * mask_ext, axis=0)
array([[4.17022005e-01, 4.17022005e-01, 4.17022005e-01],
       [8.12663088e-01, 8.12663088e-01, 8.12663088e-01],
       [1.14374817e-04, 1.14374817e-04, 1.14374817e-04],
       [1.33311280e+00, 1.33311280e+00, 1.33311280e+00]])
>>> simulation_value_ext = np.einsum('ij,i->j', e_sum_ext, weights)
array([4.36691676, 4.36691676, 4.36691676])

Here you get 3 repeated simulation values since you tiled the original rvs but if you instead tried:

>>> rvs_ext = np.random.rand(48).reshape(4,4,3)

Then you will have differ sim values.

Edit for accounting for random poisson vector.

To tweak the random lambda values try this:

>>> lambda_ext = np.array([[1,2,1,3], [2,2,3,3], [1,1,4,4]])
>>> lambda_3D = np.tile(lambda_ext.T[np.newaxis,:,:], (4,1,1)) 
>>> index_ext = np.tile(index[:,:,np.newaxis], (1,1,3))
>>> mask_3D = index_ext <= lambda_match
>>> mask_3D[:,:,0]
array([[ True,  True,  True,  True],
       [False,  True, False,  True],
       [False, False, False,  True],
       [False, False, False, False]])
>>> mask_3D[:,:,2]
array([[ True,  True,  True,  True],
       [False, False,  True,  True],
       [False, False,  True,  True],
       [False, False,  True,  True]])
>>> e_sum_ext = np.sum(rvs_ext * mask_3D, axis=0)

Now you have all skeletons to put together a script to perform this.