Birth-Death Process

4086 days ago by medlock

def steps(X, Y, *args, **kwargs): XX = [] YY = [] for i in range(len(X) - 1): XX.extend([X[i], X[i + 1]]) YY.extend([Y[i], Y[i]]) XX.append(X[-1]) YY.append(Y[-1]) return line(zip(XX, YY), *args, **kwargs) 
       
import scipy.stats def simulate(beta, mu, n0, tMin, tMax): t = tMin; N = n0; tList = [t]; NList = [N]; while t < tMax: totalRate = (beta + mu) * N dt = scipy.stats.expon(0, 1 / totalRate).rvs() if t + dt <= tMax: t += dt u = scipy.stats.uniform().rvs() if u <= beta / (beta + mu): N += 1 else: N -= 1 else: t = tMax tList.append(t) NList.append(N) return (tList, NList) 
       
beta = 0.15; mu = 0.05; N0 = 10; tMin = 0; tMax = 20; 
       
(tList, NList) = simulate(beta, mu, N0, tMin, tMax) steps(tList, NList) 
       
thMean(t) = N0 * exp((beta - mu) * t) thVar(t) = N0 * (beta + mu) / (beta - mu) * exp((beta - mu) * t) * (exp((beta - mu) * t) - 1) thStD(t) = sqrt(thVar(t)) 
       
var('t') PThStats = plot(thMean(t), (t, tMin, tMax), color = 'black') + \ plot(thMean(t) + thStD(t), (t, tMin, tMax), color = 'black', linestyle = ':') + \ plot(thMean(t) - thStD(t), (t, tMin, tMax), color = 'black', linestyle = ':') show(PThStats, axes_labels = ('$t$', '$N$')) 
       
Sims = [] PSims = [] for i in range(10): (tList, NList) = simulate(beta, mu, N0, tMin, tMax) Sims.append((tList, NList)) 
       
for i in range(len(Sims)): (tList, NList) = Sims[i] PSims.append(steps(tList, NList, hue = 0.1 * i)) show(sum(PSims) + PThStats, axes_labels = ('$t$', '$N$')) 
       
def nextTime(t, Sims): return min([min([tau for tau in tList if tau > t]) for (tList, NList) in Sims]) def getValues(t, Sims): values = [] for (tList, NList) in Sims: i = max([j for (j, tau) in enumerate(tList) if tau <= t]) values.append(NList[i]) return values 
       
t = tMin times = [t] values = getValues(t, Sims) empMeans = [mean(values)] empStds = [std(values)] while t < tMax: t = nextTime(t, Sims) times.append(t) values = getValues(t, Sims) empMeans.append(mean(values)) empStds.append(std(values)) 
       
ul = [] ll = [] for i in range(len(empMeans)): ul.append(empMeans[i] + empStds[i]) ll.append(empMeans[i] - empStds[i]) PEmpStats = steps(times, empMeans, color = 'black') + \ steps(times, ul, color = 'black', linestyle = ':') + \ steps(times, ll, color = 'black', linestyle = ':') show(sum(PSims) + PEmpStats, axes_labels = ('$t$', '$N$')) 
       
show(PEmpStats + PThStats, axes_labels = ('$t$', '$N$'))