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$'))
|
|
|
|