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

廿TT

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

窓打ち切り状況下でのポアソン過程のパラメータの最尤推定

R 確率過程

ポアソン過程

お客さんが窓口に到着する時間の間隔を確率変数 X_1,X_2,\ldots で表す.

各確率変数は独立同分布でパラメータ λ の指数分布に従うとする. このとき n \ge 1 に対して T_n = X_1 + X_2 + \cdots X_n と表記する. また T_0 = 0 とする.

 N(k) = \max \{n:T_n \le k\} はパラメータ λポアソン過程である.

状況

窓打ち切り(window censored)とは, 観測期間の幅 ( w ) が制限されている状況. 観測期間の左端点, 右端点では, 到着と到着の間隔がわからない.

f:id:abrahamcow:20150105181850p:plain

最尤推定

方針1

長さ w の窓を通してポアソン過程を観測した場合, 観測されるお客さんの人数 N(w) は平均 λwポアソン分布に従う.

\displaystyle P(N(w)=k)=\exp(- \lambda w) \frac{(\lambda w)^k}{k!}

n 回の観測があるとして, お客さんの人数の観測値を k_i ~ (i=1,\ldots,n) と表記する.

尤度は

\displaystyle L= \prod_{i=1}^{n}\exp(- \lambda w) \frac{(\lambda w)^{k_i}}{k_i!} .

対数をとると,

\displaystyle \log L = n (- \lambda w) + \sum_{i=1}^{n}\{ k_i \log (\lambda w) - \log(k_i!) \}.

λ微分して

\displaystyle \frac{d}{d \lambda}\log L = -nw + \sum_{i=1}^{n}\left\{ k_i \frac{1}{\lambda} \right\}.

0 と置いて解き, λ最尤推定量は,

\displaystyle \hat \lambda=   \frac{\sum_{i=1}^{n} k_i}{nw} .

方針2

いま, 到着と到着の間隔の時間が観測されているとする. 観測の総数を m と置く. パラメータ λ の指数分布の密度関数, 分布関数はそれぞれ f(x)=\lambda e^{-\lambda x}, F(x)=1-e^{-\lambda x} である.

到着と到着の間隔がパラメータ λ の指数分布に従う場合, 観測開始の原点 0 から次の到着までの待ち時間の分布は, やはりパラメータ λ の指数分布に従う( 再生過程における余命の分布(均衡分布) - 廿TT )ので尤度は,

\displaystyle L= \prod_{i=1}^{m} f(x_i)^{d_i} (1-F(x_i))^{1-d_i} \\
\displaystyle = \prod_{i=1}^{m} (\lambda e^{-\lambda x_i})^{d_i} (e^{-\lambda x_i})^{1-d_i}

となる. ここで d_i は右打ち切りのとき 0, 右打ち切りでないとき 1 の値を取る.

対数をとると,

\displaystyle \log L = \sum_{i=1}^{m}[ d_i  \log( \lambda e^{-\lambda x_i})+(1-d_i) \log (e^{-\lambda x}) ] \\
\displaystyle = \sum_{i=1}^{m}[ d_i ( \log \lambda  -\lambda x_i) + (1-d_i)(-\lambda x_i)]

λ微分して,

\displaystyle \frac{d}{d \lambda} \log L =\sum_{i=1}^{m}\left[ d_i \left( \frac{1}{\lambda}  -  x_i \right) + (1-d_i)(-x_i) \right] \\
\displaystyle = \sum_{i=1}^{m}\left[ \frac{d_i}{\lambda} -x_i \right]

0 と置いて解き, λ最尤推定量は,

 \displaystyle \hat \lambda= \frac{ \sum^{n}_{i=1} d_i }{ \sum^{n}_{i=1} x_i}

シミュレーション

方針 1 の推定量と, 方針 2 の推定量を比較する.

結果, どちらもほぼ同じ値になる模様.

f:id:abrahamcow:20150624055113p:plain

f:id:abrahamcow:20150624054735p:plain

シミュレーションの詳細は下記の R のコードを参照.

win=10
pics =10
simu1 <- function(){
  #set <- numeric(0)
  set <-vector("list",length(pics))
  for(j in 1:pics){
    x <- rexp(5000, 1/10) 
    t2 <- cumsum(x)
    o = runif(1,100,1000)
    sta <- which.max(t2 > o)
    en <- which.max(t2 > o+win)
    Z <- x[(sta):(en)]  
    N <- length(Z)
    A <- numeric(N)
    A[1] <- 1
    D <- rep(1,N)
    D[N] <- 0
    if(N==1){Z=win}else{
      Z[1] <- t2[sta] - o
      Z[N] <- x[en] - (t2[en] - (o+win))
    }    
    set[[j]] <-data.frame(Z,D)
  }
  return(set)
  #  return(do.call("rbind",set))
}

est1 <- est2 <-numeric(1000)

for(i in 1:1000){
  dat1 <- simu1()
  k <-sapply(dat1,function(x)dim(x)[1]) -1
  est1[i] <-sum(k)/(pics*win)
  tmp1 <- do.call("rbind",dat1)
  est2[i] <- sum(tmp1$D)/sum(tmp1$Z)
}

plot(est1,est2)
abline(0,1,col="red3",lwd=2)

boxplot(est1,est2)
abline(h=1/10, lty=2)
#MSE
> mean((est1 - 1/10)^2)
[1] 0.0009339
> mean((est2 - 1/10)^2)
[1] 0.0009339

参考文献

確率過程の基礎

確率過程の基礎