# How to simulate this Gamma expansion in a Python way

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)

• Two things: 1) Python for loops are incredibly slow; try to vectorize or use Cython instead; 2) is there any reason why you need to simulate 10 million paths..? – Helin Aug 3 '18 at 19:54
• this is part of a testing of convergence performance therefore a larger number scenarios is required – haxero Aug 11 '18 at 0:47

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)


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
>>> 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
array([[ True,  True,  True,  True],
[False,  True, False,  True],
[False, False, False,  True],
[False, False, False, False]])
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.

• how should I use index slicing? if I generate random variables like this: Nn_max = np.max(Nn) y = np.random.exponential(1,size=(nSim, Nn_max)) how could I quickly add the numbers in y according to the numbers in Nn? – haxero Aug 11 '18 at 0:53
• @haxero. amended please consider accepting the solution if it works for you – Attack68 Aug 13 '18 at 17:21
• These are very helpful. However, I have difficulty adding the third axis for each simulation. – haxero Aug 18 '18 at 8:34
• Very helpful. However, I have difficulty adding the third axis for each simulation. How to create the mask when the number_in_col is now 2-dimensional? say if the number_in_col is now: 'number_in_col = np.array([[1,2,1,3],[3,2,1,1]]) mask = number_in_col >= index ValueError: operands could not be broadcast together with shapes (2,4) (4,4) ' If I also increase the dimension of index by 1: 'index = np.tile(np.arange(1, 5)[np.newaxis,:, np.newaxis], (2,1,4)) mask = number_in_col >= index ValueError: operands could not be broadcast together with shapes (2,4) (2,4,4)' – haxero Aug 18 '18 at 8:46
• your index_ext will be 3 dim (a tiled 2 dim version), and your number_in_col will be 2 dim. I might have a look later, otherwise ask this on stackoverflow with numpy tag, someone will usually post answer quite quickly. Please accept my answer when you get it working – Attack68 Aug 18 '18 at 9:13