Python 随机模拟中的系统误差

Question

我想在python.

中使用Gillespie算法(https://en.wikipedia.org/wiki/Gillespie_algorithm)模拟一个简单的出生死亡过程

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

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，而这个模拟总是过冲——这个错误是系统性的，并且对于 [=13] 的任何选择都会保留=] 和 b。增加迭代次数不会减少此系统误差。我将其计算为

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

总是 returns 至少 10%。

我在这里错过了什么？在我的模拟中，这个系统误差的来源是什么？我是否以某种方式误解了 np.random ？我的代码一定有问题，因为错误应该缩放为 1/sqrt(# iterations)，否则。

Answer 1

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

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)

Answer 2

错误在于您检查结果的方式，特别是术语 sum(N)/len(N)。

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

from scipy.integrate import trapz

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

Python 随机模拟中的系统误差

Systematic error in Python stochastic simulation

python

numpy

stochastic

stochastic-process