廿TT

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

藤原香織、渡辺澄夫(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)

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

今回は簡単のため  \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

なんでこうなるのかは正直よくわかっていません(単にぼくがよくわかっていないという意味であり統計学的に未解決の問題ということではない)

ベイズ検定

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

変化点なしのモデルは確率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_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. アルファエラーは名目上の値と一致しました.

変化点がちょうど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のところに変化点があるのですが, 人間の目で見ても言われなければ気づかないレベルです.

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

今後の課題

普通は変化点のある位置は未知なので, ロケーションパラメータを入れて,

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

とか,

f(x_i)= a \mathrm{tanh}( x_i - b) + c

みたいな関数を考えたほうがいい気がする.

こんな風です.

y <- (Nile-mean(Nile))/sd(Nile)
x <- 1:length(y)
plot(x,y,type="l")

ll <- function(par,y,x){
  a <- par[1]
  b <- par[2]
  c <- par[3]
  fx <- a*tanh(x-c)+b
  sum(0.5*(y-fx)^2)
}
sech <- function(x)1/cosh(x)
dll <- function(par,y,x){
  a <- par[1]
  b <- par[2]
  c <- par[3]
  fx <- a*tanh(x-c)+b
  da <- -sum((y-fx)*tanh(x-c))
  db <- -sum((y-fx))
  dc <- sum((y-fx)*a*sech(x-c)^2)
  c(da,db,dc)
}
# numDeriv::grad(ll,c(2,3,4),y=y,xx=x)
# dll(c(2,3,4),y=y,x=x)

opt <- optim(c(-1,1,30),ll,dll,y=y,x=x, method = "BFGS")

plot(x,y,type="l")
lines(x,opt$par[1]*tanh(x-opt$par[3])+opt$par[2],col="royalblue", lwd=2)

f:id:abrahamcow:20200221235206p:plain

だれかやってみてはいかがですか.

その他気になる点など

検定統計量が  x_i のとり方に依存するのが若干気持ちよくない気がする.

今回はスライド(http://watanabe-www.math.dis.titech.ac.jp/users/fujiwara/doc/fujiwara_ibis2006.ppt.pdf)8ページ目に従い, -1から1の範囲に  x_i をスケーリングしたが, このスケーリング方法についてもなんらかの指針がほしい.

(これはぼくがよくわかっていないからそう感じるだけかもしれない)

アマゾン・アフィリエイトのコーナー

ベイズ統計の理論と方法

ベイズ統計の理論と方法