MATLAB代碼建模細胞複製和死亡

Question

我們從1個單元開始。它可以以速率（指數速率）1複製並以速率1死亡。設Y表示單元數。第一個事件（死亡或複製）發生在速率2.如果它是死亡 - >我們停止，因爲我們有0個單元。 如果是複製 - >我們將時間更新到t + tau，並且下一個事件現在以速率4發生（因爲2個單元格可以複製或死亡）。

由於只有2個事件可能發生，一個細胞死亡的概率是1 /（1 + 1），2個細胞的死亡概率是2 /（2 + 2）等等，複製也是一樣的。這就是爲什麼我們繪製一個從0到2的隨機數。如果這個數大於1，那麼一個細胞就會死亡，否則它會複製。

直觀地，細胞應該死時間的至少一半，因此 0細胞在時間3的概率應該是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）in平均值和std（y）/ sqrt（N）的方差。所以我會

• 在1,000次運行中返回幾次，看看結果如何變化。
• 檢查10^5運行是否表現更好。
• 準確計算std（y），或者至少將其綁定。

如果所有這些失敗，那麼它可能是一個錯誤。