廿TT

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

[Rcpp]PELT アルゴリズムで変化点の検出

今日の川柳

PELT アルゴリズム は R の changepoint パッケージですでに実装されているけど、勉強のためあらためて書いてみる。

観測値を y_{1:n}=\{y_1,\ldots,y_n\} とします。

複数の変化点があり得る場合、なんらかの情報量規準を用いて変化点検知をするには

\displaystyle G(\tau_{1:m}) = \sum_{i=1}^{m+1} C(y_{\tau_{i-1}:\tau_i}) + \beta m

を最小化する変化点 \tau_{1:m} を見つける必要があります。

ここでは、C(\cdot) は負の対数尤度とします。

\beta m はオーバーフィッティングを避けるためのペナルティ項です。

これを愚直にやると、まず 1 から n の範囲で変化点 \tau_1 を求めて、次に 1 から  \tau_1 の範囲で変化点を求めて、一方で \tau_1 から  n の範囲で変化点を求めて、……という計算を繰り返す必要があります。

PELT というアルゴリズムを使うと、n のオーダーの計算量で G(\tau_{1:m}) を最小化する \tau_{1:m} を得ることができます。

  1. F_0=-\beta, cp(0)= \rm NULL, R_1=\{0\}とする
  2. For \tau^{\ast} = 1,\ldots, n
    1. \displaystyle F_{\tau^{\ast}} = \min_{\tau < R_{\tau^{\ast}}} [ F_\tau +C(y_{(\tau+1):\tau^{\ast}})]
    2. \displaystyle \tau^{\prime} = \arg_{\tau} \{ \min_{\tau < R_{\tau^{\ast}}}  [ F_\tau +C(y_{(\tau+1):\tau^{\ast}})] \}.
    3. \displaystyle cp(\tau^{\ast})=(cp( \tau^{\prime}), \tau^{\prime})
    4. \displaystyle R_{\tau^{\ast}+1} =\{\tau \in R_{\tau^{\ast}} \cup : F(\tau) + C(y_{(\tau+1):\tau^{\ast}})\le F(\tau^{\ast})\}
  3. cp(n) に記録されている変化点を取り出す.

このアルゴリズムを Rcpp で書くとこうなりました。

#include <Rcpp.h>
using namespace Rcpp;

double cost_norm(NumericVector y, int sta, int en) {
  NumericVector suby = y[seq(sta,en)];
  double muhat = mean(suby);
  double sigmahat = sd(suby);
  return -sum(dnorm(suby,muhat,sigmahat,true));
}

// [[Rcpp::export]]
IntegerVector PELTnorm(NumericVector y){
  int n = y.length();
  NumericVector Fv(n+1);
  double pen = 2*log(n);
  Fv[0] = - pen;
  Fv[1] = 0;
  IntegerVector cpts(n);
  cpts[1] = 0;
  IntegerVector R(1);
  IntegerVector cpt_cands;
  int m;
  int tau;
  NumericVector subFv;
  LogicalVector ineq_prune;
  for(int i=1; i<n; i++){
    cpt_cands = R;
    m = cpt_cands.length();
    NumericVector seg_costs(m);
    NumericVector tmp1(m);
    for(int j=0; j<m; j++){
      seg_costs[j] = cost_norm(y, cpt_cands[j], i);
      tmp1[j] = Fv[cpt_cands[j]] + seg_costs[j] + pen;
    }
    tau = which_min(tmp1);
    Fv[i+1] = tmp1[tau];
    cpts[i] = cpt_cands[tau];
    subFv = Fv[cpt_cands];
    ineq_prune = subFv + seg_costs < Fv[i+1];
    R =  cpt_cands[ineq_prune];
    R.push_back(i);
  }
  int last = cpts[n-1];
  IntegerVector CP(1);
  CP[0]=last;
  while(last > 0){
    last = cpts[last-1];
    CP.push_back(last);
  }
  CP.sort();
  CP.erase(0);
  return CP;
}

\beta は MBIC を参考に 2\log(n) としています。

乱数でためしてみます。

library(Rcpp)
sourceCpp("PELTnorm.cpp")

plot.cps <- function(y,cps){
  plot(y,type="l")
  cps <- c(1,cps,length(y))
  means <-sapply(1:(length(cps)-1),function(i)mean(y[cps[i]:cps[i+1]]))
  for(i in 1:(length(cps)-1))segments(cps[i],means[i],cps[i+1],means[i],col="red",lwd=2)
}

set.seed(1)
y0 <-rnorm(100,0,10)
cps <- PELTnorm(y0)
plot.cps(y0,cps)

y1 <- c(rnorm(50,1,1),rnorm(50,-1,1))
cps <-PELTnorm(y1)
plot.cps(y1,cps)

y2 <- c(rnorm(25,1,1),rnorm(50,-1,1),rnorm(25,1,1))
cps <-PELTnorm(y2)
plot.cps(y2,cps)

変化点 0 のとき:
f:id:abrahamcow:20180217220421p:plain

変化点 1 つのとき:
f:id:abrahamcow:20180217220359p:plain

変化点 2 つのとき:
f:id:abrahamcow:20180217220332p:plain

参考にしたもの

PELT: Algorithm: Killick R, Fearnhead P, Eckley IA (2012) Optimal detection of changepoints with a linear computational cost, JASA 107(500), 1590–1598
MBIC: Zhang, N. R. and Siegmund, D. O. (2007) A Modified Bayes Information Criterion with Applications to the Analysis of Comparative Genomic Hybridization Data. Biometrics 63, 22-32.
GitHub - STOR-i/Changepoints.jl: A Julia package for the detection of multiple changepoints in time series.