廿TT

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

窓打ち切りデータからのワイブル再生過程のパラメータの最尤推定

再生過程

例えば電球が切れたら交換し切れたら交換し……というようなプロセスを考える.

電球それぞれの寿命 X_i>0 は独立に同分布 F(x) に従うとする. X_i は非負の連続型確率変数とする.

イベントが繰り返し生起するため, このような過程を再生過程(または更新過程; renewal process)と呼ぶ.

今回は電球の寿命にワイブル分布を仮定する.

状況

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

f:id:abrahamcow:20150105181850p:plain

最尤推定

いま, 電球の故障間隔の時間が観測されているとする. 観測の総数を n と置く. ワイブル分布の密度関数, 分布関数をそれぞれ f(x), F(x) と置く.

観測開始の原点 0 から次の故障までの待ち時間の分布は, 密度 g(x)= \frac{1-F(x)}{\mu} を持つ( 再生過程における余命の分布(均衡分布) - 廿TT ).

生存関数 G(x) = \int^{\infty}_{x} g(u) \,du についてはやや煩雑な形になるが,ガンマ関数を用いて

f:id:abrahamcow:20140303183022p:plain

と表すと数値計算上都合がいい.ここで

\displaystyle \gamma(c,x)=\int^{x}_{0}u^{c-1}e^{-u} \, du

である.

以上をすべてかけあわせて, 最大化すべき尤度は,

 \prod^{n}_{i=1} \{ f(x_i) \}^{(1-a_i)d_i} \{1-F(x_i) \}^{(1-a_i)(1-d_i)} \{g(x_i)\}^{a_i d_i} \{ G(x_i)\}^{a_i(1-d_i)}

となる. ここで d_i は右打ち切りのとき 0, 右打ち切りでないとき 1 の値を取る. a_i は原点 0 からの観測の場合に 1, そうでなければ 0 の値を取る.

シミュレーション

観測期間の幅 w を少しずつ変えて様子を見る. 赤い点線が真値.

f:id:abrahamcow:20150625003712p:plain

f:id:abrahamcow:20150625003817p:plain

f:id:abrahamcow:20150625003843p:plain

f:id:abrahamcow:20150625003910p:plain

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

pics =10
simu1 <- function(){
  #set <- numeric(0)
  set <-vector("list",length(pics))
  for(j in 1:pics){
    x <- rweibull(5000, shape=3, scale=7)
    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,A)
  }
  return(set)
  #  return(do.call("rbind",set))
}

simu1()

ll <-function(par){
#  par <- c(3,10)
  shp <- par[1]
  scl <- par[2]
  mu = scl*gamma(1+1/shp)
  tmp1 <- do.call("rbind",dat1)
  D <- as.logical(tmp1$D)
  A <- as.logical(tmp1$A)
  Z <- tmp1$Z
  lfv <- dweibull(Z[D & !A], shp, scl, log=TRUE)
  lfv = replace(lfv, is.nan(lfv), -Inf)
  lSv <- pweibull(Z[!D & !A], shp, scl, log=TRUE,lower.tail=FALSE)
  lSv = replace(lSv, is.nan(lSv), -Inf)
  lgv <- pweibull(Z[D & A],shape=shp, scale=scl, lower.tail=FALSE, log=TRUE)- log(mu)
  lgv = replace(lgv, is.nan(lgv), -Inf)
  Gv =   pgamma((Z[!D & A]/scl)^shp, shape=1+1/shp, lower.tail=FALSE)-
    (Z[!D & A]*pweibull(Z[!D & A], shape=shp, scale=scl, lower.tail=FALSE))/mu
  sum(lfv) + sum(lSv) + sum(lgv) + sum(log(Gv))
}

wins <-seq(10,100,by=10)

est_shp <- est_scl <-matrix(,100,10)

system.time({
for(j in 1:10){
  for(i in 1:100){
    win <- wins[j]
    dat1 <- simu1()
    res1 <-optim(c(1,1),ll,control=list(fnscale=-1))
    est_shp[i,j] <- res1$par[1]
    est_scl[i,j] <- res1$par[2]
  }
}})
#ユーザ   システム       経過  
#11.420      0.064     11.496 
boxplot(est_shp,xaxt="n",xlab="window size", main="shape parameter")
axis(side=1, labels=wins, at=1:10)
abline(h=3, lty=2, lwd=2, col="red3")

boxplot(est_scl,xaxt="n",xlab="window size", main="scale parameter")
axis(side=1, labels=wins, at=1:10)
abline(h=7, lty=2, lwd=2, col="red3")

plot(wins,apply((est_shp-3)^2,2,mean),type="o",
     xlab="window size",ylab="MSE", main="shape parameter")

plot(wins,apply((est_scl-7)^2,2,mean),type="o",
     xlab="window size",ylab="MSE",main="scale parameter")