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

廿TT

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

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

R 仮説検定 確率過程

モデル

変化点のあるポアソン分布についてはいろんな研究がなされていますが、一番単純と思われる手法を試します。

解析対称はイギリスの炭鉱事故の発生件数のデータです。

3. Tutorial — PyMC 2.3.6 documentation から取得しました。

#R のコード
dat <- c(4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
          3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
          2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
          1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
          0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
          3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
          0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1)
plot(dat,type="h",ylab="frequency")

f:id:abrahamcow:20151029003418p:plain

40年ごろを境に事故の発生件数が減少していることが伺えるので、その変化を抽出できるようなモデルを考えます。

事故の発生件数はポアソン分布に従うと仮定し、ある時期を境にポアソン分布の平均が変化したと考えると良さそうです。

そこで i 年目の事故の発生件数を x_i とし、x_i はパラメータ \theta_iポアソン分布に従うというモデルを立てました(i=1,\ldots,T)。

ただし \theta _tt <\tau のとき \theta_1、そうでなければ \theta_2 の値を取ります。

このモデルのもとで、尤度関数は下のようになります。

\displaystyle L_1 = \prod_{i=1}^{\tau} \frac{\theta _1^{x_i}e^{-\theta_1}}{x_i !} \prod_{j=\tau+1}^{n} \frac{\theta _1^{x_i}e^{-\theta_1}}{x_i !}

対数をとり整理すると、

\displaystyle \log L_1 = S_{\tau} \log \theta _1 +(S_T-S_\tau) \log \theta _2 - \theta _1 \tau - \theta _2 (T-\tau) -\sum_{i=1}^{T} \log (x_i !)

となります。
ただしここで S_t=x_1 + \cdots + x_t です。

\log L_1 \theta _1\theta_2最尤推定量をそれぞれ代入して、

\displaystyle l_1(\tau) = S_T \log(S_T/\tau)\\
\displaystyle ~~~~+(S_T-S_\tau)\log\{ (S_T-S_\tau)/(T-\tau\}-S_T-\sum_{i=1}^{T}\log (x_i !)

 l_1\tau について最大化することで、変化点の推定値が求まります。

コーディング

which.max 関数を使い、総当りで  l_1 を最大にする \tau を求めることにします。

x <- dat
n <- length(x)
ST <- sum(x)
l1_tau <- function(tau){
  Stau <- sum(x[1:tau])
  Stau*log(Stau/tau) + (ST-Stau)*log((ST-Stau)/(n-tau))
}
tau <-which.max(sapply(1:(n-1),l1_tau))

\tau は41でした。

> tau
[1] 41

 \theta _1\theta_2最尤推定値はそれぞれ、3.01、0.91でした。

> theta1 <-mean(x[1:tau])
> theta2 <-mean(x[(tau+1):n])
> theta1
[1] 3.097561
> theta2
[1] 0.9142857

f:id:abrahamcow:20151128152033p:plain

plot(x,type="h",lwd=2)
curve(ifelse(x>tau,theta2,theta1),add=TRUE,col="royalblue",lwd=2)

尤度比検定

変化点がないと仮定した場合の対数尤度、

\displaystyle l_0 = S_T \log (S_T/T) -S_T - \sum_{i=1}^{T} \log( x_i!)

l_1 の差を評価することで、尤度比検定を行うことができます。

-2 (l_0 -l_1) は漸近的に自由度 3 - 1 = 2 のカイ二乗分布に従うと言われています。

3 - 1 = 2 は変化点のあるモデルのパラメータ数 - 変化点のないモデルのパラメータ数です。

l0 <-ST*log(ST/n) - ST - sum(log(factorial(x)))
Stau <- sum(x[1:tau])
l1 <- Stau*log(Stau/tau) + (ST-Stau)*log((ST-Stau)/(n-tau)) - ST - sum(lfactorial(x))

p 値は以下のように求まります。

> pchisq(-2*(l0-l1),2,lower.tail = FALSE)
[1] 1.418789e-15

有意水準を 5% と決めていたとすると、「パラメータを増やすことによって尤度が改善していない」という帰無仮説が棄却され、変化点があるとみなしたほうがいいことがわかります。

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

参考文献

入門・演習 数理統計

入門・演習 数理統計