Python随机模拟中的系统误差

Question

我想用python中的Gillespie算法（https://en.wikipedia.org/wiki/Gillespie_algorithm）模拟一个简单的生灭过程。你知道吗

在每一瞬间，每个人都有出生的概率和死亡的概率。我相信下面的代码应该模拟种群的动态：

import numpy as np 
a = 10.0 # birth rate
b = 1.0 # death rate
n = 0.0 # initial population
t = 0.0 # initial time
T = [] # times of births or deaths
N = [] # population timeseries
for _ in range(1000000): # do 1000000 iterations
    P = [a, b*n] # rate of birth and death
    P0 = sum(P) # total rate
    # step the time
    t = t + (1/P0)*np.log(1/np.random.random())
    # select the transition
    R = P[0]/P0 # probability of birth
    r = np.random.random() # choose a random number
    # enact the transition
    if r<R: # birth
        n = n+1
    elif r>=R: # death 
        n = n-1
    N.append(n) # update the output lists
    T.append(t)

这具有预期的行为：绘制N与T显示具有明确统计模式的随机总体。然而，我有一个严重的困惑。你知道吗

这个模型的解析解说N的平均值应该是a/b，而这个模拟总是超调——这个误差是系统性的，并且为a和b的任何选择所保留。增加迭代次数并不能减少这种系统误差。我计算为

(sum(N)/len(N)-a/b)/(a/b)*100 # error in the mean value in percent

总回报至少10%。你知道吗

我错过了什么？在我的模拟中，这个系统误差的来源是什么？我是不是误解了np.random？我的代码一定有问题，因为错误应该扩展为1/sqrt（#迭代），否则。你知道吗

Answer 1

错误在于如何检查结果，特别是术语sum(N)/len(N)。你知道吗

您必须随着时间的推移进行整合：

from scipy.integrate import trapz

theo = a/b
obs = trapz(N,T)/T[-1]
(obs-theo)/theo
# 0.0029365091018275892

Answer 2

你需要考虑时间吗？我不太熟悉这些东西，但这个简化版本至少给出了随着时间的推移接近于零的结果。你知道吗

import random

TIME = 100000.0

a = 10.0  # birth rate
b = 1.0   # death rate
n = 0     # initial population
t = 0.0   # initial time

t_average_pop = 0.0

while True:
    P0 = a + b * n

    # step the time
    step = random.expovariate(P0)
    t += step
    t_average_pop += n * step

    if t >= TIME:
        break

    # select the transition
    if random.uniform(0, P0) < a:
        # birth
        n += 1
    else:
        # death
        n -= 1

average_pop = t_average_pop / t
expected = a / b

print((average_pop - expected) / expected * 100)