MATLAB 代码建模细胞复制和死亡

Question

我们从 1 个单元格开始。它可以以速率（指数速率）1 复制并以速率 1 死亡。让 Y 表示细胞数。第一个事件（死亡或复制）以速率 2 发生。如果它是死亡 -> 我们停止，因为我们有 0 个细胞。 如果是复制 -> 我们将时间更新为 t+tau，下一个事件现在以速率 4 发生。（因为 2 个细胞可以复制或死亡）。

由于只能发生 2 个事件，因此一个细胞的死亡概率为 1/(1+1)，2 个细胞为 2/(2+2)，以此类推，复制也是如此。这就是为什么我们从 0 到 2 中抽取一个随机数。如果这个数字 >1，那么一个细胞就会死亡，否则它会复制。

直观地说，至少有一半的时间细胞应该死亡，因此 在时间 3 出现 0 个细胞的概率应该是 P(Y=0)>0.5（实际上 答案是 3/4）。但是当我将此代码放入 for 循环并运行它时 1000 次，我得到 Y=0 的次数约为 400，即 0.4

t=0;
rr=1; %rate of replication 1 cell -> 2 cells
rd=1; %rate of death 1 cell -> 0 cells
Y=1; %initial number of cells
while t<3 && Y>0 % interested in probabilities of number of cells at time t=0,t=1,t=2,t=3
    r = 2*rand; %draws a random number from 0 to 2
    tau=exprnd(2*Y); %since the total rate of all possible events is replication+death=2 for each cell
    if t+tau < 3 %if the event happens before 3 seconds
        if r>1 %death
        Y=Y-1;
        else Y=Y+1; %otherwise replication
        end 
    elseif t+tau > 3 %if the next event happens after 3 seconds, we are not interested.
        Y; 
        t; 
        break
        end;
    t=t+tau; %update time from t to t+tau
end

Answer 1

好吧，在运行调试过程之前，请检查您的统计信息！

根据中心极限定理，如果我们将 y 定义为随机变量，当 Y=0 时 y=1，否则 y=0，则 y 的 N 次运行的平均值应收敛于均值 P(Y=0) , 和 std(y)/sqrt(N) 的方差。所以我会

• 在 1000 次运行中返回几次，看看结果如何变化。
• 检查 10^5 次运行是否表现更好。
• 精确计算 std(y)，或者至少绑定它。

如果所有这些都失败了，那么它可能是一个错误。