藤原香織、渡辺澄夫（2006）ベイズ尤度比検定による変化点検出のシミュレーション

これです：

http://watanabe-www.math.dis.titech.ac.jp/users/fujiwara/doc/fujiwara_ibis2006.ppt.pdf

この研究成果は論文化されているようですが, オープンアクセスではないみたいです：

https://search.ieice.org/bin/summary.php?id=j91-d_4_889&category=D&year=2008&lang=J&abst=

本文中のコードはすべてRです.

カイ二乗検定

データのある点 $x_i$ は固定とし $y_i$ が正規分布,

$y_i \sim \mathcal{N}(f(x_i), \sigma^ 2)$

$f(x_i)= a \mathrm{tanh}(b x_i)$

に従うというモデルを考えます.

つまり変化点のある系列データに, こんな感じの曲線を当てはめるようなモデルを考えています.

f:id:abrahamcow:20200221235206p:plain

今回は簡単のため $\sigma^ 2 = 1$ になるよう, データが適切にスケーリングされているとします.

変化点なしのモデルは $a=0$ , $b=0$ に相当します.

ふつうの尤度比検定を行うことを考えます. 対数尤度比は, 変化点ありのモデルの尤度を変化点なしのモデルの尤度でわり算し,

$l_n = -\sum_i ( ( \hat f(x_i) )^ 2-2y_i \hat f(x_i))/2$

$\hat f(x_i)= \hat a \mathrm{tanh}(\hat b x_i)$

です. $\hat a$ , $\hat b$ は最尤推定量です.

このモデルでは, 帰無仮説のもとでの対数尤度比の2倍がカイ二乗分布に従いません.

シミュレーションしてみます.

sech <- function(x){
  1/cosh(x)
}
ll <- function(par,x,y){
  a <- par[1]
  b <- par[2]
  fx <- a*tanh(b*x)
  sum(fx^2-2*y*fx)/2
}
dll <- function(par,x,y){
  a <- par[1]
  b <- par[2]
  fx <- a*tanh(b*x)
  da <- sum((fx-y)*tanh(b*x))
  db <- sum((fx-y)*x*a*sech(b*x)^2)
  c(da,db)
}

n <- 100
x <- seq(-1,1,length.out = n)
stat <- numeric(10000)
for (i in 1:10000) {
  y <- rnorm(n)
  opt <- optim(runif(2,-1,1),ll,dll,x=x,y=y,method = "BFGS")
  stat[i] <- -2*opt$value
}
hist(stat,breaks = "FD",freq = FALSE)
curve(dchisq(x,df=1),add=TRUE,n=1000,col="red")
pv <- pchisq(stat,df=1,lower.tail = FALSE)
hist(pv)

検定統計量のヒストグラムを書いてみると一見カイ二乗分布っぽい.

f:id:abrahamcow:20200221224957p:plain

でも p 値を出してみると小さい値が出やすいことがわかる.

f:id:abrahamcow:20200221225047p:plain

print(round(mean(pv<0.05),2))
# 0.07

なんでこうなるのかというと, パラメータの a, b でそれぞれ符号を入れ替えても得られる関数 $f(x_i)$ の形はまったく一緒で, 最尤推定量がユニークじゃないからみたいです（正直まだよく理解できてないので, そのうち整理します.）

ベイズ尤度比検定

対立仮説と帰無仮説を確率分布で設定します.

変化点なしのモデルは確率1で $a=0$ , $b=0$ が出る分布（デルタ分布）として, 変化点ありのモデルは $a$ と $b$ がともに区間 [-1, 1] の一様分布から生成されることを考えます.

パラメータ $a$ と $b$ をまとめて $w$ とおくと, ベイズ尤度比は

$L_n(Y) = \int \exp( -\sum_i ( ( f(x_i) )^ 2-2y_i f(x_i))/2 ) / 4 \,dw$

です. $1/4$ は一様分布の密度です.

$f(x)$ を2次の項まで0の周りでテイラー展開すると,

$\frac{\partial}{\partial a} f(x) = \mathrm{tanh}(b x)$

$\frac{\partial}{\partial b} f(x) = a \mathrm{sech}^ 2(b x) x$

$\frac{\partial ^2}{\partial a ^2} f(x) = 0$

$\frac{\partial ^2}{\partial b ^2} f(x) = a 2 \mathrm{sech}(b x) (- \mathrm{tanh}(b x)) x ^2$

$\frac{\partial ^2}{\partial a \partial b} f(x) = \mathrm{sech}^ 2(b x) x$

より

$f(x) \approx abx$

と近似できますので,

$L_n(Y) \approx \int \exp( -\sum_i ( \frac{1}{2}( a b x_i )^ 2 - a b x_i y_i) ) / 4 \,dw$

が成り立ちます.

$N = \frac{1}{2n}\sum_i x_i ^ 2$ ,

$M = \frac{1}{\sqrt{n}} \sum_i x_i y_i$

とおくと, $N$ は定数で確率変数ではありません. $M$ は（帰無仮説のもとで）平均 0, 分散 $2N$ の正規分布に従います.

$L(M) = \int \exp( -( a b ) ^2 N + a b \sqrt{n}M) ) / 4 \,dw$

とおき, デルタ関数（ $\delta(x)$ ）を使って変形すると,

$L(M) = \int _ {-1} ^ {1} \int _ {-1} ^ {1} \int _ {-1} ^ {1} \delta(t-ab) \exp( -t ^ 2 N + t \sqrt{n}M) ) / 4 \,dt \, da \, db$

です. $t'= \sqrt{n}t$ とおくと, 積分範囲は $[-1, 1]$ から $[-\sqrt{n}, \sqrt{n}]$ に変わり,

$L(M) = \int _ {-\sqrt{n}} ^ {\sqrt{n}} \int _ {-1} ^ {1} \int _ {-1} ^ {1} \delta(t/\sqrt{n}-ab) \exp( -t ^ 2 N + t M) ) / 4 \,dt \, da \, db$

となります.

ここで, （なぜ成り立つかはよくわかってないのですが）

$\int _ {-1} ^ {1} \int _ {-1} ^ {1} \delta(t/\sqrt{n}-ab) \, da \, db = -2 \log(|t|/\sqrt{n})$

をそういうものだと思って使い,

$L(M) = \int _ {-\sqrt{n}} ^ {\sqrt{n}} \\{ -2 \log(|t|/\sqrt{n}) \\} \exp( -t ^2 N + t M) ) / 4 \,dt$

と変形すると,

$L(M) = \int _ {0} ^ {\sqrt{n}} \\{ - \log(|t|/\sqrt{n}) \\} \exp( -t ^2 N) \\\ \times \\{ \exp( t M)+\exp( -(t M) )) \\} /2 \,dt \\\ = \int_0 ^ {\sqrt{n}} \frac{1}{\sqrt{n}} (-\log (|t|/\sqrt{n}) )e ^ {-Nt ^2} \mathrm{cosh} (Mt) \, dt$

という式が出てきて, $M$ が正規分布なので, 有意水準 5% とすると約標準偏差2つ分の範囲の外側が棄却域になります.

つまり検定統計量が

$L(2\sqrt{2 N}) = \int_0^{\sqrt{n}} \frac{1}{\sqrt{n}} (-\log (|t|/\sqrt{n}) ) e ^{-Nt ^2} \mathrm{cosh} (2\sqrt{2 N}t) \, dt$

を超えたら棄却します.

アルファエラーのシミュレーション

シミュレーションしてみます.

積分は数値積分すると遅いので（ぼくは割とせっかちなので）モンテカルロ積分しています.

まずは帰無仮説のもとでのシミュレーション.

library(parallel)
L <- function(par,x,y){
  a <- par[1]
  b <- par[2]
  fx <- a*tanh(b*x)
  exp(-sum(fx^2-2*y*fx)/2)
}

# library(cubature)
# intL <- adaptIntegrate(L, lowerLimit = c(-1, -1), upperLimit = c(1, 1),
#                        x=x,y=y)
# intL$integral/4

x <- seq(-1,1,length.out = n)
simstat <- function(i,beta=0,x){
  set.seed(i)
  n <- length(x)
  y <- c(rnorm(n/2,-beta/2),rnorm(n/2,beta/2))
  wsamp <- matrix(runif(2*1000,-1,1),1000,2)
  return(-log(mean(apply(wsamp, 1, L,x=x,y=y))))
}
out_alpha <- mclapply(1:10000,simstat,x=x,
                      mc.cores = detectCores())

N <- mean(x^2)/2
z_a <- qnorm(0.975)
integrand <- function(t,N,n,z_a){
  (-log(abs(t)/sqrt(n)))*exp(-N*t^2)*cosh(z_a*sqrt(2*N)*t)
}
int <- integrate(integrand,0,sqrt(n),N=N,n=n,z_a=z_a)
thres <- log(n)/2 - log(int$value)

hist(unlist(out_alpha))
abline(v=thres)

これが検定統計量の分布です.

f:id:abrahamcow:20200221233102p:plain

round(mean(unlist(out_alpha)<thres),2)
# 0.05

棄却される確率は約 0.05.
アルファエラーは名目上の値と一致しました.

n を減らしてやってみます.

n = 10 のとき.

n <- 10
x <- seq(-1,1,length.out = n)
simstat <- function(i,beta=0,x){
  set.seed(i)
  n <- length(x)
  y <- c(rnorm(n/2,-beta/2),rnorm(n/2,beta/2))
  wsamp <- matrix(runif(2*1000,-1,1),1000,2)
  return(-log(mean(apply(wsamp, 1, L,x=x,y=y))))
}
out_alpha <- mclapply(1:10000,simstat,x=x,
                      mc.cores = detectCores())

N <- mean(x^2)/2
integrand <- function(t,N,n,z_a){
  (-log(abs(t)/sqrt(n)))*exp(-N*t^2)*cosh(z_a*sqrt(2*N)*t)
}
int5 <- integrate(integrand,0,sqrt(n),N=N,n=n,z_a=qnorm(0.975))
thres5 <- log(n)/2 - log(int5$value)
int10 <- integrate(integrand,0,sqrt(n),N=N,n=n,z_a=qnorm(0.95))
thres10 <- log(n)/2 - log(int10$value)
int1 <- integrate(integrand,0,sqrt(n),N=N,n=n,z_a=qnorm(0.995))
thres1 <- log(n)/2 - log(int1$value)

round(mean(unlist(out_alpha)<thres5),2)
round(mean(unlist(out_alpha)<thres10),2)
round(mean(unlist(out_alpha)<thres1),2)

有意水準	アルファエラー
5%	4%
10%	9%
1%	1%

サンプルサイズが10もあれば通例使われてる有意水準で, アルファエラーが名目上の数字と一致しています.

検出力のシミュレーション

変化点がちょうど0のところにあり, 平均が -0.5/2 から 0.5/2 に変化するとき,

out_05 <- mclapply(1:10000,simstat,x=x,beta=0.5,
                      mc.cores = detectCores())

hist(unlist(out_05))
abline(v=thres)

検定統計量の分布はこれです.

f:id:abrahamcow:20200221233501p:plain

round(mean(unlist(out_05)<thres),2)
# 0.59

検出力は約0.59でした.

変化点がちょうど0のところにあり, 平均が -1/2 から 1/2 に変化するとき,

out_1 <- mclapply(1:10000,simstat,x=x,beta=1,
                   mc.cores = detectCores())

hist(unlist(out_1))
abline(v=thres)

検定統計量の分布はこれです.

f:id:abrahamcow:20200221233713p:plain

round(mean(unlist(out_1)<thres),2)
# 0.99

検出力は約0.99でした.

標準偏差1で変化の量が1の時系列というのは例えばこんな感じです.

#unlist(out_1)[1]<thres
#[1] TRUE
set.seed(1)
y <- c(rnorm(50,-1/2),rnorm(50,1/2))
plot(x,y,type="l")

f:id:abrahamcow:20200221233902p:plain

0のところに変化点があるのですが, 人間の目で見ても言われなければ気づかないレベルです.

検出力はめちゃくちゃ高いといえるでしょう.

変化点のロケーションがずれている場合

対立仮説のモデルは変化点の位置が固定されていたので, そこの指定が間違っていたときにどれくらい検出力が下がるのかを調べます.

平均を -1/2 から 1/2 に変化させ, 時系列の長さを100とし, 変化点のロケーションを系列データの 10個目, 20個目, 30個目, 40個目, 50個目とずらしてみます.

f:id:abrahamcow:20200222155454p:plain

変化点が端っこにある場合はやはりちょっと難しいようですが, それでもけっこう検出できてますね.

R のコードを貼ります.

library(parallel)
L <- function(par,x,y){
  a <- par[1]
  b <- par[2]
  fx <- a*tanh(b*x)
  exp(-sum(fx^2-2*y*fx)/2)
}
n <- 100
x <- seq(-1,1,length.out = n)
simstat <- function(i,beta=0,x,cp){
  set.seed(i)
  n <- length(x)
  y <- c(rnorm(cp,-beta/2),rnorm(n-cp,beta/2))
  wsamp <- matrix(runif(2*1000,-1,1),1000,2)
  return(log(mean(apply(wsamp, 1, L,x=x,y=y))))
}
N <- mean(x^2)/2
z_a <- qnorm(0.975)
integrand <- function(t,N,n,z_a){
  (-log(abs(t)/sqrt(n)))*exp(-N*t^2)*cosh(z_a*sqrt(2*N)*t)
}
int <- integrate(integrand,0,sqrt(n),N=N,n=n,z_a=z_a)
thres <- -log(n)/2 + log(int$value)

out_10 <- mclapply(1:10000,simstat,x=x,beta=1,cp=10,
                  mc.cores = detectCores())
out_20 <- mclapply(1:10000,simstat,x=x,beta=1,cp=20,
                  mc.cores = detectCores())
out_30 <- mclapply(1:10000,simstat,x=x,beta=1,cp=30,
                   mc.cores = detectCores())
out_40 <- mclapply(1:10000,simstat,x=x,beta=1,cp=40,
                   mc.cores = detectCores())
out_50 <- mclapply(1:10000,simstat,x=x,beta=1,cp=50,
                   mc.cores = detectCores())

beta <- numeric(5)
beta[1] <- mean(unlist(out_10)>thres)
beta[2] <- mean(unlist(out_20)>thres)
beta[3] <- mean(unlist(out_30)>thres)
beta[4] <- mean(unlist(out_40)>thres)
beta[5] <- mean(unlist(out_50)>thres)
plot(1:5*10,beta,type = "b",
     xlab = "location of change point",
     ylab = "power")

その他気になる点など

検定統計量は $x_i$ のとり方に依存するのかな？

今回はスライド（http://watanabe-www.math.dis.titech.ac.jp/users/fujiwara/doc/fujiwara_ibis2006.ppt.pdf）8ページ目に従い, -1から1の範囲に $x_i$ をスケーリングした.

もっと学びたい人のための論文リンク

特異な場合という問題意識が欠如してる可能性もありそうですね．

混合正規分布の均一性検定も特異モデルの検定ですが，ベイズ検定と変分ベイズ検定が研究されてますねhttps://t.co/2NevGQD5pT https://t.co/gylfxs1S18

さっきの変化点検知と共通の著者がいることからも界隈は広くないです
— カイヤン🦈博論執筆Chap5 (@389jan) 2020年7月16日