読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

出生死滅過程(birth-death process)のシミュレーション

出生死滅過程の性質

※これは定義ではない.

パラメータ λ, μ の出生死滅過程 {X(t),t \ge 0} において, 次の出生または死亡までの時間は, パラメータ λ + μ の指数分布に従う.

状態 i から i-1 へ推移する確率は μ/( λ + μ ), i+1 へ推移する確率は λ/( λ + μ ) である.

R によるシミュレーション

Ninit <- 10  ##Initial population size
mu <- 0.3   ##death rate
lambda <- 0.7 ##birth rate
BirthDeathSim <-function(){
    Ts <- rexp(1000,rate=mu+lambda)
    CTs <- cumsum(Ts)
    Ts <- Ts[CTs < 10]
    CTs <- CTs[CTs < 10]
    n <-length(CTs)
    tmp <-sample(c(1,-1),n,prob=c(lambda/(mu+lambda),mu/(mu+lambda)), replace=TRUE)
    length(CTs)
    CTs <- c(0,CTs)
    path <-cumsum(c(Ninit,tmp))
    data.frame(time=CTs,populations=path)
}

set.seed(8)
dat <-BirthDeathSim()

library(ggplot2)
theme_set(theme_bw(15))
ggplot(dat) +
  geom_step(aes(time,populations))

f:id:abrahamcow:20150815174453p:plain

図は乱数で生成したサンプルパス.

推定

n をジャンプの総数, B を出生イベントの総数, D を死亡イベントの総数として,  \hat \lambda = B/n,  \hat \mu = D/n最尤推定量のような気がするが, よくわからない.

Brate <- Drate <- numeric(10000)

for(i in 1:10000){
  dat <-BirthDeathSim()
  BorD <-diff(dat$populations)
  n <- length(BorD)
  Drate[i]<-sum(BorD == -1)/n
  Brate[i]<-sum(BorD == 1)/n
}

ggplot(data.frame(value=Drate)) +
  geom_histogram(aes(x=value),position="identity")+
  geom_vline(xintercept =0.3,colour="red") +
  labs(y=expression(hat(mu)))

ggplot(data.frame(value=Brate)) +
  geom_histogram(aes(x=value),position="identity")+
  geom_vline(xintercept =0.7,colour="red") +
  labs(y=expression(hat(lambda)))

f:id:abrahamcow:20150815175137p:plain

f:id:abrahamcow:20150815175148p:plain

図は推定値の分布を示したヒストグラム. 赤い線は真値.

> summary(Drate)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.2000  0.2857  0.2967  0.4000  1.0000       2 
> summary(Brate)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.6000  0.7143  0.7033  0.8000  1.0000       2 

参考文献

確率モデル入門

確率モデル入門