Rで解析：時系列の変化点を検出する「scanr」パッケージ

長期間にわたって記録した時系列データでは、ある時点を境に平均や分布そのものが変わることがあり、その位置を正確に把握することが分析の出発点になります。しかし、変化した位置を目視で見極めたり、窓幅を変えながら検定を繰り返したりする作業には手間がかかります。

「scanr」パッケージは、ノンパラメトリックな検定を用いて単変量の時系列から変化点を検出するパッケージです。複数の窓幅による走査結果を投票でまとめる関数、単一の窓幅で走査する関数、CUSUM統計量やワッサースタイン統計量で変化位置を求める関数が収録されています。また、検出した変化点と真の変化点を突き合わせて精度指標を算出することや、時系列と変化点を重ねて可視化することも可能です。

パッケージバージョンは0.1.1。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~

データ分析の一連の流れを手を動かして学ぶ定番書。ただしコードはRではなくPython中心なので、その点は注意。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。作図に用いるggplot2などは依存パッケージとして自動的に導入されます。

# パッケージのインストール
install.packages("scanr")

# パッケージの読み込み
library("scanr")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

変化点とは、時系列の途中で統計的な性質（平均・分散・分布の形）が切り替わる位置を指します。scanrは、時系列の各位置について、その前後の区間を比べる統計量を計算し、値が大きく跳ね上がる位置を変化点の候補とします。区間の比較には、平均のずれに反応するCUSUM統計量と、分布の形の違いに反応するワッサースタイン統計量が用意されています。比較に使う区間の長さ（窓幅）は結果に影響するため、複数の窓幅で走査し、候補を投票でまとめる関数も用意されています。

以降の例では、京都府和束町の茶園に設置した地温センサーが、120日間にわたって記録した日ごとの平均地温（摂氏）を題材にします。61日目に畝を覆う被覆資材を切り替えたため、そこを境に地温の水準がおよそ3度上がった、という状況を想定した擬似データを次のコードで作成します。

# 乱数の種を固定して結果を再現できるようにする
set.seed(2026)

# 被覆資材の切り替え前、60日分の日平均地温
before <- rnorm(60, mean = 15, sd = 1.2)

# 切り替え後、60日分の日平均地温
after <- rnorm(60, mean = 18, sd = 1.2)

# 前後をつないで120日分の時系列にする
soil_temp <- c(before, after)

[スキャン窓サイズの自動選択]：default_window_sizesコマンド

系列の長さから、複数窓での走査に使う窓幅の並びを自動的に決めます。下限と上限のあいだから、指定した個数の窓幅を等間隔に選びます。scan_cpdにwindow_sizesを渡さない場合は、内部でこの関数と同じ規則が使われます。

# 120点の時系列に対する既定の窓幅の並び
default_window_sizes(120L)
[1] 15 18 19 20 22


# 窓幅の下限・上限と個数、乱数の種を指定して求める
default_window_sizes(
  n = 120L,
  min_window = 12L,
  max_window = 48L,
  n_windows = 6L,
  seed = 2024L
)
[1] 13 22 28 40 43 45

[単変量時系列における変化点の検出]：scan_cpdコマンド

複数の窓幅で時系列を走査し、それぞれの結果を投票でまとめて変化点を返します。本パッケージの中心となる関数です。change_typeで、平均のずれ・分散のずれ・分布全体の違いのいずれに反応させるかを選びます。

# 3種類の窓幅でアンサンブルの走査を実行する
fit <- scan_cpd(
  soil_temp,
  window_sizes = c(20L, 30L, 40L),
  random_state = 2026L,
  change_type = "mean"
)

# 検出結果の概要を表示する
fit
scanr change-point result
observations: 120
change points: 60


# 変化点の位置だけを取り出す
fit$change_points
[1] 60

被覆資材を切り替えた61日目の直前、60番目の観測が変化点として検出されました。

＜おすすめのRに関する書籍です＞

RStudioではじめるRプログラミング入門 | ギャレット・グロ-ルマンド, 大橋真也, 長尾高弘

統計の前に、まずRをプログラミング言語として覚える一冊。身近な題材で関数・条件分岐・ループを学べる。

[単一の窓幅に対する走査の実行]：scan_single_windowコマンド

窓幅を1つだけ指定して走査します。投票をはさまないため、その窓幅での検定結果がそのまま得られます。窓幅による結果の変わり方を確かめたいときに使います。

# 窓幅を30に固定して単一窓の走査を実行する
sw <- scan_single_window(
  soil_temp,
  window_size = 30L,
  random_state = 2026L,
  change_type = "mean"
)

# 検出結果を表示する
sw
scanr single-window result
window size: 30
change points: 60

[CUSUM統計量による平均変化の特定]：ts_cusumコマンド

累積和（CUSUM）にもとづく統計量が最大になる位置を返します。平均が一段階ずれる形の変化に向いた、単一の走査です。引数は数値ベクターのみです。

# CUSUM統計量で平均が変化した位置を求める
ts_cusum(soil_temp)
[1] 60

[ワッサースタイン統計量による分布変化の特定]：ts_wassersteinコマンド

前後の区間の分布間のワッサースタイン距離にもとづく統計量から、変化点の位置と各位置の統計量を返します。平均だけでなく、ばらつきや分布の形の違いにも反応します。返り値はリストで、change_pointに変化点の位置が入ります。

# ワッサースタイン統計量で分布が変化した位置を求める
tw <- ts_wasserstein(soil_temp)
# 変化点の位置を取り出す
tw$change_point
[1] 60

＜おすすめのRに関する書籍です＞

超入門! Rでできるビジュアル統計学 学会・論文発表に役立つデータ可視化マニュアル

学会発表や論文の図表づくりに寄せた実践書。伝わるグラフをRで作るための手順をまとめている。

www.amazon.co.jp

[1次元ワッサースタイン距離の計算]：one_wasserstein_distanceコマンド

2つの数値ベクターを標本とみなし、その経験分布のあいだの1次元ワッサースタイン距離（1-Wasserstein距離）を返します。ts_wassersteinが内部で用いる、区間どうしの距離そのものを取り出す関数です。

# 切り替え前後それぞれ30日分の地温を取り出す
before30 <- soil_temp[1:30]
after30 <- soil_temp[91:120]

# 2つの期間の分布間の1次元ワッサースタイン距離を求める
one_wasserstein_distance(before30, after30)
[1] 2.812299

[変化点検出の精度評価]：cpd_metricsコマンド

真の変化点と推定した変化点を突き合わせ、適合率・再現率・F1値と、区間の重なりにもとづくcovering指標をまとめて返します。toleranceで、何点までのずれを一致とみなすかを決めます。ここでは、別の手法が31番目と60番目を候補として返したと仮定し、scan_cpdの検出位置を真の変化点として評価します。

# 別の手法が返した変化点の候補（誤検出を1つ含む）
candidate <- c(31L, 60L)

# scan_cpdの検出位置を真の変化点とみなして候補を評価する
cpd_metrics(
  true_cps = fit$change_points,
  estimated_cps = candidate,
  n = 120L,
  tolerance = 5L
)
$matches
  true estimated distance
1   60        60        0

$precision
[1] 0.5

$recall
[1] 1

$f1
[1] 0.6666667

$covering
[1] 0.7583333

60番目は許容範囲内で一致したため再現率は1、いっぽう31番目は誤検出のため適合率は0.5となっています。

[時系列データと変化点の可視化]：vis_time_seriesコマンド

時系列を折れ線で描き、指定した変化点の位置に縦線を重ねたggplotオブジェクトを返します。true_change_pointsに真の変化点を渡すと、検出位置と並べて確認できます。x_label・y_labelで軸ラベルを指定できます。

# 地温の推移に、検出した変化点を重ねて描画する
vis_time_series(
  soil_temp,
  change_points = fit$change_points,
  x_label = "経過日数",
  y_label = "日平均地温"
)

単一の窓幅で走査したときと複数の窓幅で走査したときとで結果を見比べたい場合は、それぞれのchange_pointsをこの関数に渡して並べて描くと、窓幅の選び方による違いを確認できます。

＜おすすめのRに関する書籍です＞

RとPythonで学ぶ[実践的]データサイエンス&機械学習【増補改訂版】 | 有賀友紀, 大橋俊介

RとPythonを併走させながらデータサイエンスと機械学習を学ぶ。両方の言語を行き来する人に向く。

この記事が誰かの役に立ちますように。