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

廿TT

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

ワイブル過程のパラメータの最尤推定

ワイブル過程とは

ワイブル過程とは強度(intensity)関数に \lambda(t)=\beta m t^{m-1} を仮定した非定常ポアソンのことで, おもに信頼度成長のモデルに使われる.

この強度関数はワイブル分布のハザード関数と一致する.

ワイブル分布とは別ものなので紛らわしい.

power-law ポアソン過程とか, ワイブル−ポアソン過程とか呼ばれることもあるらしい.

  • m >1 のとき, 事象の起こりやすさ(強度)は時間とともに増加する.
  • m=1 のとき, 事象の起こりやすさは時間に依存しない. 点の集合は定常ポアソン過程になる. 点と点の間隔は指数分布に従う.
  • m<1 のとき, 事象の起こりやすさは時間とともに減少する.

最尤推定

時刻 t における強度が  \lambda(t) の非定常ポアソン過程を, 区間 (0, T] で観測した場合の尤度関数は以下のようになる.

\displaystyle L(\beta,m)= \prod_{i=1}^{n} \lambda ( t_{i}) \exp\left(-\Lambda(T)\right)

ただし,
\displaystyle \Lambda(T) = \int_{0}^{T}\lambda\left( u\right) \,du.

ワイブル過程の場合, パラメータ m,\beta最尤推定量は閉じた形で求まる.

\displaystyle \log L(\beta,m)= \sum \log\lambda ( t_{i}) -\Lambda(T) \\
=\displaystyle \sum_{i=1}^{n}\left( \log \beta + \log m + (m-1)\log(t_i) \right) - \beta T^m \\
=\displaystyle n\log \beta + n\log m + (m-1)\sum_{i=1}^{n}\log(t_i) - \beta T^m

これを微分して 0 と置いて解く.

β については,

\displaystyle \frac{\partial}{\partial \beta} \log L(\beta,m)= \frac{n}{\beta} -  T^m = 0,

\displaystyle \hat \beta=n/T^m.

m については,

\displaystyle \frac{\partial}{\partial m} \log L(\beta,m) = \frac{n}{m}+\sum_{i=1}^{n}\log t_i -\beta T^m \log T=0

 \hat \beta を代入して,

\displaystyle \frac{n}{m}+\sum_{i=1}^{n}\log t_i -n \log T=0,

\displaystyle \hat m=\frac{n}{\sum_{i=1}^{n}\log (t_i/T_0)}.

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

ワイブル過程の乱数を生成して, その乱数からパラメータを推定することで, 推定量の精度を見る.

乱数の生成に関しては, RK's Musings: Simulating NHPP in RGenerating a non-homogeneous Poisson process | Freakonometrics を参考にした.

愚直なコードなのでけっこう遅い.

rWP <- function(beta,m,t_max){
  t <- 0
  s <- 0
  v <- seq(0,t_max, by = 0.01)
  Lambda <- function(t)beta*(t^m)
  Lambda_inv <- function(s){
    min(v[Lambda(v)>=s]) #Λ(t)>=tを満たす最小のt
  }
  X <- numeric(0)
  while(t <= t_max){
    u <- runif(1)
    s <- s -log(u)
    t <- Lambda_inv(s)
    X <- c(X, t)
  }
  X[-length(X)]
}
dat <-rWP(1.5,0.5,1000)

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

res1 <- matrix(nrow=100,ncol=2)
colnames(res1) <- c("beta","m")
for(i in 1:100){
  dat <- rWP(2,2,10)
  res1[i,] <- estimator(dat,10)
}

res2 <- matrix(nrow=100,ncol=2)
colnames(res2) <- c("beta","m")
for(i in 1:100){
  dat <- rWP(1.5,0.5,1000)
  res2[i,] <- estimator(dat,1000)
}

\beta =2,m=2,T=10 を仮定したとき

> summary(res1)
      beta              m        
 Min.   :0.8717   Min.   :1.697  
 1st Qu.:1.4031   1st Qu.:1.971  
 Median :1.8369   Median :2.033  
 Mean   :1.8796   Mean   :2.049  
 3rd Qu.:2.1627   3rd Qu.:2.150  
 Max.   :4.3186   Max.   :2.316 

\beta =1.5,m=0.5,T=1000 を仮定したとき

> summary(res2)
      beta              m         
 Min.   :0.2666   Min.   :0.3496  
 1st Qu.:0.8919   1st Qu.:0.4648  
 Median :1.3055   Median :0.5145  
 Mean   :1.4529   Mean   :0.5245  
 3rd Qu.:1.9311   3rd Qu.:0.5717  
 Max.   :4.2900   Max.   :0.7566