Figure 10.15:

Probability density for single particle stochastic growth model at different times.

Code for Figure 10.15

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# Converted from singleparticle.m - Poisson distribution for single particle
import numpy as np
from scipy.special import gammaln
from misc import octave_save

tfin = 1.0
nts  = 5
tvec = np.linspace(0.1, 0.9, nts)

# large lambda
lam_big = 10000
nsim_big = int(tfin * lam_big)
igridbig = np.arange(1, nsim_big + 1)
pbig = np.zeros((nsim_big, nts))
for i, tau_frac in enumerate(tvec):
    tau   = tfin * tau_frac
    logp  = -lam_big*tau + (igridbig - 1)*np.log(lam_big*tau) - gammaln(igridbig)
    pbig[:, i] = np.exp(logp)

# small lambda
lam_small = 100
nsim_small = int(tfin * lam_small)
igridsmall = np.arange(1, nsim_small + 1)
psmall = np.zeros((nsim_small, nts))
for i, tau_frac in enumerate(tvec):
    tau   = tfin * tau_frac
    logp  = -lam_small*tau + (igridsmall - 1)*np.log(lam_small*tau) - gammaln(igridsmall)
    psmall[:, i] = np.exp(logp)

table1 = np.column_stack([igridsmall, psmall])
table2 = np.column_stack([igridbig,   pbig])
octave_save('singleparticle.dat', ('table1', table1), ('table2', table2))