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

廿TT

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

未知の変化点があるモデルでは AIC が使えない

R 時系列

モデル

時系列データ x_t (t=1,\ldots, n) があるとします.
このデータが, 変化点( \tau)以前では平均 \mu _1, 標準偏差 1 の正規分布に従い, 変化点から後には平均 \mu_2, 標準偏差 1 の正規分布に従うと考えます.
標準偏差は既知とします.

 \displaystyle x_t =\begin{cases} \mu_1 + \varepsilon _t & t \le \tau \\ \mu_2 + \varepsilon _t & t > \tau \end{cases}

ここで \varepsilon _t は標準正規分布に従う確率変数です.

変化点 \tau最尤推定するには, 対数尤度関数に \tau が与えられたときの \mu _1, \mu _2最尤推定量(標本平均)を代入して尤度が最大になる点を探してやればよさそうです.

最大化すべき対数尤度関数は以下です.

 \displaystyle l(\tau)= \sum_{i=1}^{\tau} (\log(\phi (x_i | \hat \mu_1))+ \sum_{j=\tau+1}^{n} (\log(\phi (x_j | \hat \mu_2)).

ここで \phi標準偏差 1 の正規分布の密度関数,
 \displaystyle \hat \mu_1 = \frac{1}{\tau}\sum_{i =1}^{\tau} x_i,  \displaystyle \hat \mu_2 = \frac{1}{n-\tau}\sum_{i = \tau+1}^{n} x_i
です.

R で推定

乱数で適当なデータを作って, \tau, \mu _1, \mu _2 を推定してみます.

\tau=50, \mu _1=-1, \mu _2=1, n=100と設定しました.

set.seed(1)
x=c(rnorm(50,-1),rnorm(50,1)) #データの生成
ll1_f <- function(tau,n,x){ #尤度関数の定義
  sum(dnorm(x[1:tau],mean(x[1:tau]),log=TRUE))+
    sum(dnorm(x[(tau+1):n],mean(x[(tau+1):n]),log=TRUE))
}
l1v <- sapply(1:99,ll1_f,n=100,x=x) #尤度が最大になる点を総当りで探す
cp <-which.max(l1v)
mu1hat = mean(x[1:cp])
mu2hat = mean(x[(cp+1):100]) 
plot(x,type="b")
segments(x0=0,y0=mu1hat,x1=cp,y1=mu1hat,col="red",lwd=2)
segments(x0=cp+1,y0=mu2hat,x1=100,y1=mu2hat,col="red",lwd=2)

f:id:abrahamcow:20170214010307p:plain

無事, 変化点と \mu _1, \mu _2 が推定されました.

次に AIC で変化点ありのモデルと変化点なしのモデルを比較してみます.

変化点なしのモデルは,

 \displaystyle x_t = \mu_0 + \varepsilon _t

です.

変化点なしのモデルのパラメータは  \mu_0 の 1 個, 変化点ありのモデルのパラメータは \tau, \mu _1, \mu _2 の 3 個です.

ll1 <-max(l1v) #変化点ありのモデルの対数尤度
ll0 <- sum(dnorm(x,mean(x),log=TRUE)) #変化点なしのモデルの対数尤度
-2*ll0+2 #AIC
-2*ll1+2*3

変化点なしのモデルの AIC は367.345, 変化点ありのモデルの AIC は 269.65 となり, 変化点ありのモデルのほうが AIC が小さい=良いモデルとして選択されました。

ここからが本題です。

シミュレーション

あえて変化点なしのモデルで乱数を生成して, 変化点ありのモデルとなしのモデルで AIC を比較してみます.

真のモデルは変化点なしです.

どちらの AIC が小さくなるでしょうか.

simAIC <- function(i,n){
  x=rnorm(n)
  l1v <- sapply(1:(n-1),ll1_f,n=n,x=x)
  cp <-which.max(l1v)
  ll1 <-max(l1v)
  ll0 <- sum(dnorm(x,mean(x),log=TRUE))
 c(AIC0 = -2*ll0 + 2,
   AIC1 = -2*ll1 + 2*3)
}

library(parallel)
n10 <-simplify2array(mclapply(1:10000,simAIC,n=10,mc.cores=detectCores()-1))
n100 <-simplify2array(mclapply(1:10000,simAIC,n=100,mc.cores=detectCores()-1))
n1000 <-simplify2array(mclapply(1:10000,simAIC,n=1000,mc.cores=detectCores()-1))

mean(n10["AIC0",] < n10["AIC1",])
mean(n100["AIC0",] < n100["AIC1",])
mean(n1000["AIC0",] < n1000["AIC1",])

サンプルサイズ n を 10, 100, 1000 と変化させて, それぞれ 1 万回のシミュレーションをおこないました.

変化点なしのモデルが選択された割合は以下のようになりました.

 n 割合
10 76%
100 47%
1000 28%

サンプルサイズが増えても, 正しいモデルが選択されない, というかサンプルサイズが増えるほど, 変化点ありのモデルが選ばれやすくなることがわかります.

なんでこうなるのか, ではどうしたらいいのかはそのうち書きます.