廿TT

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

変化点のあるポアソン過程のパラメータの最尤推定

尤度関数

変化点のあるポアソン分布のパラメータの最尤推定 - 廿TT では、生起したイベントの個数に着目しましたが、生起の間隔に着目してモデル化することもできます。

 \Delta t_i = t_i - t_{i-1} とすると、変化点のない強度(intensity)λ の定常ポアソン過程では、点と点の間隔は指数分布に従うので、対数尤度を、

 \displaystyle L(\Delta t_1, \ldots ,\Delta t_n; \lambda) \\
 \displaystyle =\log \prod_{i=2}^{n} \lambda e^{-\lambda \Delta t_i}\\
 \displaystyle =(n-1)\log \lambda - \lambda (t_n-t_1)

と書くことができます。

最尤推定\hat \lambda = (n-1)/(t_n-t_1) を代入して、

 \displaystyle L_0 = (n-1)\log \frac{n-1}{t_n-t_1}- (n-1)

となり、これが最大化された対数尤度です。

さて、ある時期を境にポアソン過程のレートパラメータが変化するとします。

変化点がひとつのとき、対数尤度は、

 \displaystyle L_1 = L(\Delta t_1, \ldots ,\Delta t_{n_1}; \lambda _1) \\
 \displaystyle ~~+ L(\Delta t_{n_1 + 1}, \ldots ,\Delta t_n; \lambda _2) \\
\displaystyle = (n_1-1) \log \frac{n_1-1}{t_{n_1}-t_1} \\
\displaystyle ~~+ (n-n_1) \log \frac{n-n_1}{t_{n}-t_{n_1}} - (n-1)

となり、これを最大化する  n_1 を求めることで変化点が求まります。

同様に、変化点がふたつのとき、対数尤度は、

 \displaystyle L_2 = L(\Delta t_1, \ldots ,\Delta t_{n_1}; \lambda _1)\\
\displaystyle ~~+L(\Delta t_{n_1+1}, \ldots ,\Delta t_{n_2}; \lambda _2) \\
\displaystyle ~~+ L(\Delta t_{n_2 + 1}, \ldots ,\Delta t_n; \lambda _3) \\
\displaystyle = (n_1-1) \log \frac{n_1-1}{t_{n_1}-t_1} \\
\displaystyle~~+ (n_2-n_1) \log \frac{n_2-n_1}{t_{n_2}-t_{n_1}} \\
\displaystyle~~+(n-n_2) \log \frac{n-n_2}{t_{n}-t_{n_2}} - (n-1)

となります。

R を用いた例題

イギリスの炭鉱事故のデータを使います。boot パッケージの coal データには1851年から1962年までの事故の発生日が記録されています。

f:id:abrahamcow:20151229233939p:plain

data("coal",package = "boot")
ti <-unlist(coal)
n <- length(ti)
plot(ti,1:n,type="s")

l1_tau <- function(n1){
  (n1-1)*log((n1-1)/(ti[n1]-ti[1])) + 
    (n-n1)*log((n-n1)/(ti[n]-ti[n1])) - (n-1)
}
l1v <-sapply(1:n,l1_tau)
l1 <- max(l1v,na.rm = TRUE)
n1 <-which.max(l1v)

l0 <- (n-1)*log((n-1)/(ti[n]-ti[1]))-(n-1) 

pchisq(-2*(l0-l1),2,lower.tail = FALSE)
#3.426829e-16 
####
l2_tau <- function(n1,n2){
  (n1-1)*log((n1-1)/(ti[n1]-ti[1])) + 
    (n2-n1)*log((n2-n1)/(ti[n2]-ti[n1])) +
    (n-n2)*log((n-n2)/(ti[n]-ti[n2]))- (n-1)
}

n <- length(ti)
mat1 <-matrix(,n-1,n-1)

for(j in 2:(n-1)){
  for(i in 2:(n-1)){
    if(i<j) mat1[i,j] <- l2_tau(i,j)
  }
}
mat1 <-ifelse(mat1==Inf,NA,mat1)
l2 <-max(mat1,na.rm = TRUE)
n12 <-which(mat1==l2,arr.ind = TRUE)

plot(ti,1:n,type="s")
abline(v=ti[n12],col="royalblue",lty=2,lwd=2)

尤度比検定の結果、変化点がひとつのモデルは 5% 水準で棄却され、変化点はふたつあるとみなしたほうがよさそうです。

> pchisq(-2*(l1-l2),2,lower.tail = FALSE)
[1] 0.005070965

ただし上記の推定法では帰無分布は自由度1のカイ二乗分布の最大値の分布になるため、この検定は正確ではありません。

f:id:abrahamcow:20151229234447p:plain

定常ポアソン過程の累積ハザード関数 \Lambda(t)=\lambda t を重ねてプロットしてみます。

n1 <- n12[1]
n2 <- n12[2]
lambda1 <- (n1-1)/(ti[n1]-ti[1])
lambda2 <- (n2-n1)/(ti[n2]-ti[n1])
lambda3 <- (n-n2)/(ti[n]-ti[n2])
lambda <- function(t){
  if(t<ti[n1]){
    lambda1*(t-ti[1]) + 1
  }else if(t<ti[n2]){
    lambda2*(t-ti[n1])+n1
  }else{
    lambda3*(t-ti[n2])+n2
  }
}
lambda <-Vectorize(lambda)
plot(ti,1:n)
curve(lambda(x),add=TRUE,col="red",lwd=3)

f:id:abrahamcow:20151230072007p:plain

各期間の λ を表にまとめておきます。

tab1 <-data.frame(c(ti[1],
ti[n1],
ti[n2]),
c(lambda1,
lambda2,
lambda3))
colnames(tab1) <-c("年","平均発生回数")
rownames(tab1) <-NULL

λ は単位時間(ここでは年)あたりのイベントの発生回数と解釈できます。

平均発生回数
1851.20〜1890.19 3.18
1890.19〜1947.69 1.08
1947.69〜 0.28