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

廿TT

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

観測期間にギャップのあるワイブル過程のパラメータの最尤推定

R 確率過程

最尤推定

f:id:abrahamcow:20151108220440p:plain

ワイブル過程のパラメータの最尤推定 - 廿TT の続きです.

Crow (1988) に出てくる例題をやります.

x_iを区間 (0, S_1] でのイベントの発生時刻, y_iを区間 [S_2, T] でのイベントの発生時刻として,非定常ポアソン過程を区間 (0, S_1][S_2,T]S_1 \le S_2)で観測した場合の尤度関数は以下のようになります.

\displaystyle L(\beta,m)= \prod_{i=1}^{n_1} \lambda ( x_{i}) \exp\left(-\Lambda(S_1)\right) \\
~~~~\prod_{j=1}^{n_2} \lambda (y_{j})\exp\left( -( \Lambda(T)-\Lambda(S_2) ) \right)

対数をとり,

\displaystyle \log L(\beta,m)= \sum_{i=1}^{n_1} \log\lambda (x_{i}) -\Lambda(S_1) \\
~~~~+\sum_{j=1}^{n_2} \lambda (y_{j}) -( \Lambda(T)-\Lambda(S_2) )

ワイブル過程の強度関数 \lambda(t)=\beta m t^{m-1} を代入して,

\log L(\beta,m)\\
 = \sum_{i=1}^{n_1}( \log \beta+ (m-1) \log x_i )  -\beta S_1^m  \\
~~~~+ \sum_{j=1}^{n_2} (\log \beta+ (m-1) y_{j}) -( \beta T^m -\beta S_2^m )\\
 = n_1 \log \beta + (m-1)\sum_{i=1}^{n_1}\log x_i\\
~~~~  -\beta S_1^m + n_2 \log \beta + (m-1)\sum_{j=1}^{n_2} y_{j} -\beta(  T^m - S_2^m )

これを微分して 0 と置いて解きます. 閉じた形では推定量が求まりません.

β については,

\displaystyle \frac{\partial}{\partial \beta} \log L(\beta,m)= \frac{n_1}{\beta} -S_{1}^{m}+\frac{n_2}{\beta} -  (T^m -S^{m}_{2})= 0,

\displaystyle \hat \beta=\frac{n_1+n_2}{T^m - S^{m}_{2}}.

m については,

\displaystyle \frac{\partial}{\partial m} \log L(\beta,m) = \\
\displaystyle  \frac{n_1}{m}+\sum_{i=1}^{n}\log x_i -\beta S_{1}^{m} \log S_1+\frac{n_2}{m}\\
 +\sum_{j=1}^{n_2}\log y_j -\beta( T^{m} \log T-S_{2}^{m} \log S_{2}^{m})  =0,

\displaystyle \hat m=(n_1+n_2) \big/ \left( \hat \beta (S_{1}^{\hat m} \log S_1 + T^{\hat m}\log T- S^{\hat m}_{2} \log S_2) \\
~~~~ +\sum_{i=1}^{n_1}\log x_i + \sum_{j=1}^{n_2}\log y_j \right) .

 \hat \lambda\hat m を更新式としてイテレーションを回す必要があります.

R による計算例

使用するデータはこれです.

ex2 <- c(.5, .6, 10.7, 15.6, 18.3, 19.2, 19.5, 25.3, 39.2, 39.4, 43.2, 44.8, 47.4, 65.7,
 88.1,  97.2, 104.9, 105.1, 120.8, 195.6, 217.1, 219.0, 257.5, 260.4, 281.3,
283.7, 289.8, 306.6, 328.6, 357.0, 371.7, 374.7, 393.2, 403.2, 466.5, 500.9,
501.5,  518.4, 520.7, 522.7, 524.6, 526.9, 527.8, 533.6, 536.5, 542.4,543.2,
545.0,  547.4, 554.0, 554.1, 554.2, 554.8, 556.5, 570.6, 571.4, 574.9, 576.8,
578.8, 583.4, 584.9, 590.6, 596.1, 599.1, 600.1, 602.5, 613.9, 616.0, 616.2,
617.1, 621.4, 622.6, 624.7, 628.8,  642.4, 684.8, 731.9, 735.1,753.6, 792.5,
803.7, 805.4, 832.5, 836.2, 873.2, 975.1)

(600,625) の期間は欠測しているとします.

S1=500; S2=625; t_max <- 1000
x <-ex2[ex2<=S1]
y <-ex2[ex2>=S2]
estimator2 <- function(x,y,beta0=1,m0=1,maxit=100){
  maxit <- 1000
  n <- length(x)+length(y)
  for(i in 1:maxit){
    mhat <- n/(beta0*((S1^m0)*log(S1)+(t_max^m0)*log(t_max)-(S2^m0)*log(S2))-
                 (sum(log(x))+sum(log(y))))
    betahat <- n/(S1^mhat+t_max^mhat-S2^mhat)
    if(abs(mhat-m0)<1e-6 & abs(betahat-beta0)<1e-6)break
    m0<-mhat
    beta0<-betahat
  }
  c("beta"=betahat,"m"=mhat)
}

論文に書いてある推定値は \lambda = 1.108, m= 0.559.

> estimator2(x,y)
     beta         m 
1.1081554 0.5592317 

計算結果が一致しました.

完全データを使って求めた場合とだいぶ値違うけどいいのかこれ.

> estimator <- function(dat){
+   n <-length(dat)
+   mhat <- n/sum(log(t_max/dat))
+   betahat <- n/(t_max^mhat)
+   c("beta"=betahat,"m"=mhat)
+ }
> estimator(ex2)
     beta         m 
0.4534404 0.7593261 

f:id:abrahamcow:20151110120629p:plain

青い線が完全データから求めた場合, 赤い線が欠測データから求めた場合, 黒い点線が欠測期間.

fit_comp <- estimator(ex2)
fit_gap <-estimator2(x,y)
Lambda <- function(t,beta,m)beta*(t^m)
plot(ex2,1:length(ex2),xlab="time",ylab="cumulative number of events")
curve(Lambda(x,fit_comp[1],fit_comp[2]),add=TRUE,col="royalblue",lwd=2)
curve(Lambda(x,fit_gap[1],fit_gap[2]),add=TRUE,col="tomato",lwd=2)
abline(v=S1,lty=2)
abline(v=S2,lty=2)

おまけ

Example 1 のデータもタイプしてしまったので置いておきます.

ex1 <- c(.5, .6, 10.7, 15.6, 18.3, 19.2, 19.5, 25.3, 39.2, 39.4, 43.2, 44.8, 47.4, 65.7,
88.1, 97.2, 104.9, 105.1, 120.8, 195.6, 217.1, 219.0, 257.5, 260.4, 281.3,
283.7, 289.8, 306.6, 328.6, 357.0, 371.7, 374.1, 393.2, 403.2, 466.5, 520.7,
556.5, 571.4, 621.4, 628.8, 642.4, 684.8, 731.9, 735.1, 753.6, 792.5, 803.7,
805.4, 832.5, 836.2, 873.2, 975.1)

論文に書いてある推定値は \lambda = 1.132, m= 0.554.