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

Question

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

由于只有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=1和Y=0分别定义y为带有y=0的随机变量时，则Y的N次平均运行应收敛到平均P(Y=0)，而std(y)/sqrt(N)则为方差。所以我会

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

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