RでChangeFinder - J's blog

変化点検出のアルゴリズムとして、ChangeFinderという有名なアルゴリズムがあります。Rによる実装例が見つからなかったので実装しました。今回は１次元のデータに対する変化点検出を考えます。アルゴリズムの詳細は記述しませんので、後述の参考群を参照してください。

実装

まず、準備として予め２つの関数を定義しておきます。
１つ目は、アルゴリズム内で用いる平滑化用の関数です。
２つ目は、擬似逆行列を計算する関数です。計算上、逆行列が存在しない場合もありますので、擬似逆行列を計算する必要があります。今回はムーア・ペンローズ逆行列を実装しています。

## 平滑化用関数
MW <- function(dat, width, na.rm = TRUE) {
  N <- length(dat)
  x <- numeric(N)
  for(i in 1L:(width-1L)) {
    x[i] <- mean(dat[1L:i], na.rm = na.rm)
  }
  for(i in width:N) {
    x[i] <- mean(dat[(i-width+1L):i], na.rm = na.rm)
  }
  return(x)
}
## 擬似逆行列計算用関数（MASSパッケージより）
ginv <- function (X, tol = sqrt(.Machine$double.eps)) {
  if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X))) 
    stop("'X' must be a numeric or complex matrix")
  if (!is.matrix(X)) 
    X <- as.matrix(X)
  Xsvd <- svd(X)
  if (is.complex(X)) 
    Xsvd$u <- Conj(Xsvd$u)
  Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
  if (all(Positive)) 
    Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
  else if (!any(Positive)) 
    array(0, dim(X)[2L:1L])
  else Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))
}

ChangeFinderの要の部分になります。ARで定めた時系列モデルからの逸脱度を計算します。さらにこのアルゴリズムはオンラインアルゴリズムであるため、忘却係数を用いたパラメータの更新も行います。スコア計算には対数損失を用いています。

SDAR1 <- function(x, ar_order, forgetting, mu = 0, co = numeric(ar_order+1), w = numeric(ar_order), S = 0) {
  N <- length(x)
  score <- numeric(N)
  for(t in (ar_order+1L):N) {
    mu <- (1-forgetting)*mu + forgetting*x[t]
    co <- (1-forgetting)*co + forgetting*(x[t]-mu)*(x[t-(0:ar_order)]-mu)
    w <- c(ginv(matrix(co[abs(rep(1:ar_order, ar_order)-rep(1:ar_order, each=ar_order))+1], ar_order, ar_order)) %*% co[-1])
    x_hat <- c(w %*% (x[t-(1:ar_order)]-mu)) + mu
    S <- (1-forgetting)*S + forgetting*(x[t]-x_hat)^2
    score[t] <- (log(2*pi*S) + (x[t] - x_hat)^2 / S) / 2
  }
  list(Score = score, x = x)
}

いよいよChangeFinderアルゴリズムの完成になります。SDARと平滑化を繰り返すのみですので、ここは簡潔に済みます。
またここで、パラメータの第1Step・第2Stepの初期値をどうするかという問題があるのですが今回はすべて0としています。

CF <- function(x, ar_order, forgetting, width1, width2 = width1, na.rm = TRUE) {
  Step1 <- SDAR1(x, ar_order, forgetting)$Score
  Step2 <- MW(Step1, width1, na.rm)
  Step3 <- MW(SDAR1(Step2, ar_order, forgetting)$Score, width2, na.rm)
  return(Step3)
}

適用例

シミュレーションにはF# でオンライン変化点検知 - 捨てられたブログの例を利用しています。

x <- c(rnorm(300, 0.7, 0.05), 
       rnorm(300, 1.5, 0.05), 
       rnorm(300, 0.6, 0.05), 
       rnorm(300, 1.3, 0.05))
CF_test <- CF(x, 1, 0.01, 7, 7)
plot(x, col = rep(c(1:2, 1:2), each = 300), ylim = range(x, CF_test))
lines(CF_test, col = 4)