2025-04-21

時系列データにおける多重構造変化検出

はじめに

時系列データ解析では，「ある時点で回帰関係が突然変化し，レジームが切り替わる」現象が多く観察されます。たとえばマクロ経済の景気転換点や金融市場のボラティリティ変動，政策介入後の効果検証などが典型例です。本記事では，R の strucchange パッケージに含まれる breakpoints() 関数による、複数箇所の構造変化点を同時に推定する枠組みを解説します。

本記事は以下の論文を参考文献としています。

https://onlinelibrary.wiley.com/doi/10.1002/jae.659

詳しいパッケージについては以下を参照してください。

cran.r-project.org

モデルと変数の定義

観測時刻を $t=1,2,\dots,n$ とし，説明変数ベクトルを $x_t\in\mathbb{R}^p$ ，応答変数を $y_t$ ，誤差項を $u_t$ とします。
変化点を $T_1,\dots,T_m$

ただし、以下の制約を満たすものとします。

$1 \le T_1 \lt \cdots \lt T_m \lt n$

と仮定すると，データは $m+1$ のセグメントに分割され，各セグメント $j$ （ $T_{j-1}+1$ から $T_j$ ）では別個の回帰係数 $\beta_j$ に従うモデルとなります。すなわち

$y_t = x_t^\top \beta_j + u_t,\quad t = T_{j-1}+1,\dots,T_j,\quad j = 1,\dots,m+1$

ただし $T_0=0,;T_{m+1}=n$ と定義します。

最小二乗法による最適化

未知量として各セグメントの回帰係数 $\beta_1,\dots,\beta_{m+1}$ と変化点位置 $T_1,\dots,T_m$ を同時に推定するため，各時刻の残差平方和（RSS）を

$\mathrm{RSS} = \sum_{j=1}^{m+1}\sum_{t=T_{j-1}+1}^{T_j} \bigl(y_t - x_t^\top \beta_j\bigr)^2$

と定義し，これを最小化するパラメータを

$\left(\hat\beta_1,\dots,\hat\beta_{m+1},\hat T_1,\dots,\hat T_m\right)$
$= \operatorname*{arg\,min}_{\substack{\beta_1,\dots,\beta_{m+1}\\1\le T_1<\cdots<T_m<n}}$
$\mathrm{RSS}$

として求めます。

全探索による組み合わせ数と計算量

変化点位置を全組み合わせ探索するには，仕切り位置候補が $n-1$ 箇所あり，そこから $m$ 箇所を選ぶ組み合わせ数

$\binom{n-1}{m}$

をすべて評価する必要があります。さらに変化点数 $m$ を固定せずに $m=0,1,\dots,n-1$ をすべて試す場合には，組み合わせ数の合計が

$\sum_{m=0}^{n-1}\binom{n-1}{m} = 2^{n-1}$

となり，計算量は指数関数オーダー

$O\bigl(2^{n-1}\bigr)$

に膨れ上がります。観測点 $n$ が数百を超えるケースでは，全探索は困難です。

解決策：動的計画法

以上のように，全組み合わせ探索では計算量が $O(2^{n-1})$ に爆発し，実用上の限界を迎えます。そこで Bai and Perron（2003）が提案した動的計画法（Dynamic Programming）を用いることで，前もってすべての区間 $[i,j$ ;(1\le i<j\le n)] の残差和を $\tfrac{n(n-1)}{2}$ 件計算・蓄積し，Bellman の原理に基づく再帰的結合で最適解を多項式時間

$O(n^2)$

で得ることが可能になります。詳細なアルゴリズムは本記事の範囲外としますが，これにより指数爆発する全探索と異なり，実務的に十分高速な推定が実現できます。

動的計画法について理解が深まれば，別途記事にまとめます。

2025-03-12

Hyndman–KhandakarアルゴリズムによるARIMAモデルの自動選定

はじめに

本記事では、時系列予測において広く利用されるARIMAモデルの自動選定手法として、Hyndman–Khandakarアルゴリズムについて解説します。

この記事は以下を参考文献としています。

www.jstatsoft.org

ARIMAモデルのパラメータ（自己回帰成分 $p$ 、差分次数 $d$ 、移動平均成分 $q$ など）は基本的に、推定を行う前に選定する必要があります。

しかし、Hyndman–Khandakarアルゴリズムでは、統計的検定と情報量基準を組み合わせることで、データに応じた最適なモデルを自動的に決定できます。

この手法は、Rの $auto.arima()$ 関数などに実装され、専門知識がなくても高精度な予測が可能となっています。　

ARIMAモデルについて

ARIMAモデルは、時系列データのパターンを解析し将来の値を予測するための代表的な手法です。基本的な非季節性のARIMAモデルは、自己回帰（AR）、差分（I）、移動平均（MA）の3つの主要な成分で構成され、次のように表現されます。

$\phi(B)(1-B)^d y_t = c + \theta(B)\epsilon_t$

ここで、 $y_t$ は観測された時系列データを示し、 $B$ はバックシフト演算子です。 $\phi(B)$ は自己回帰成分を表す多項式、 $\theta(B)$ は移動平均成分を示す多項式、 $d$ は系列の差分次数、 $c$ は定数項、そして $\epsilon_t$ はホワイトノイズとして扱われる誤差項を意味します。

ARIMAモデルでは自己回帰成分の次数 $p$ 、移動平均成分の次数 $q$ 、および差分次数 $d$ を適切に選定する必要があります。

このパラメータの選定を人の手で行う場合、数多くの可能性の中から最適な組み合わせを見つけ出す必要があり、容易ではありません。

Hyndman–Khandakarアルゴリズムは、統計的検定と情報量基準を組み合わせることで、最適なパラメータの組み合わせを自動的に選定できる手法として提案されました。

Hyndman–Khandakarアルゴリズムの詳細解説

Hyndman–Khandakarアルゴリズムは、まず系列の定常性を評価し、必要な差分回数を自動的に決定するところから始まります。

非季節性の時系列の場合、元の系列に対して単位根検定であるKPSS検定(ちなみにKPSSは考案者4人の頭文字だそう)を行います。

ここで帰無仮説が棄却された場合、系列が非定常であると判断され、1階差分 $\Delta y_t = y_t - y_{t-1}$ を適用します。

この検定を、系列が定常とみなされるまで繰り返すことで、最小限の差分次数 $d$ が自動的に決定されます。

※ここでKPSS検定は必ずしも単位根が存在することを意味するわけでは無いことに注意されたい

次に、差分後の定常系列に対して、自己回帰成分と移動平均成分の次数を決定する段階に入ります。すべての可能なモデルの組み合わせを評価するのは現実的ではないため、アルゴリズムはまず初期モデルを設定します。

たとえば、非季節データでは $ARIMA(2, d, 2)$ を初期候補として採用します。初期モデルから出発し、パラメータ $p, q$ を個々に ±1 ずつ変更した近傍モデル群を生成し、各モデルの尤度に基づく情報量基準（ $AIC$ )で評価します。もし近傍モデルの中でAICが改善されるものがあれば、そのモデルが新たな基準モデルとして採用され、再びその近傍モデルを探索するプロセスが繰り返されます。こうして、探索はAICがこれ以上改善されなくなるまで続けられ、最終的に最適なパラメータの組み合わせが自動的に選択されます。

※以下自信なし

さらに、モデルの推定過程では、ARおよびMA多項式の根の大きさが評価され、因果性や逆行性を損なわないような制約が課せらるそうです(あまり追えていません)。

具体的には、多項式の根の絶対値が1に近すぎる場合、そのモデルは不安定であるとみなされ候補から除外されるそう。また、非線形最適化の過程でエラーが発生したモデルも採用されないため、安定して推定可能なモデルのみが最終候補として残る仕組みになるらしいです、、、

こうした工夫により、Hyndman–Khandakarアルゴリズムは、差分次数の自動決定と局所的なパラメータ探索を組み合わせることで、手作業での試行錯誤を極力排除し、堅牢かつ効率的に最適なARIMAモデルを構築することを可能にしているそうです、なんかわからんけど凄い？

少し考えたらわかる気もしなくもないですが、一旦置いておきます。またしっかり追えたら追記します。

最後にMermaidで作成したグラフを添付します。

この一連の自動モデル選定と予測の流れにより、ユーザーは専門的なパラメータ選定の知識がなくても、信頼性の高い時系列予測を実現できるようになります。

R の $auto.arima()$ 関数にて実装されているので多くの人が知らずに用いていると思います。

ありがとう、技術。

以上、今回はHyndman–KhandakarアルゴリズムによるARIMAモデルの自動選定についてまとめました。

2025-03-08

Interrupted Time Series Analysisにおける多重共線性の考察

Interrupted Time Series Analysisでよく利用される segmented regression のパラメータについて、なぜ従来のモデル（Equation 1）が多重共線性の問題を引き起こすのか、そしてその問題をどのように解決するのかを数学的な観点から解説します。

以下のアンサー記事のようなものになっています。

qiita.com

本記事作成にあたって以下の論文を参考にしています。(この論文には多重共線性ではなく、モデルのバイアスについて記載されています)

pubmed.ncbi.nlm.nih.gov

1. モデルの定式化と問題の概要

まず、一般的なモデル設定として、次の Equation 1 を考えます。

$Y_t = \beta_0 + \beta_1 T + \beta_2 X_t + \beta_3 (T \times X_t)$

ここで、

$T$ は経過時間を表します。
$X_t$ は介入前は 0、介入後は 1 となるダミー変数です（介入開始時刻を $T_i$ とする）。
$T \times X_t$ は $T$ と $X_t$ の交互作用項です。

2. 多重共線性とランク落ちの数学的意味

2.1 交互作用項の性質

介入前は $X_t = 0$ のため、交互作用項は

$T \times X_t = 0$

となりますが、介入後 $X_t = 1$ では

$T \times X_t = T \times 1 = T$

となります。
すなわち、介入後のデータにおいては、変数 $T$ と交互作用項は常に同じ値となります。
より形式的には、介入後の観測に対して $Z = T \times X_t$ とおくと、

$Z = T$

となります。このため、モデル内でこれら二つの変数は同じ情報を持つことになり、数学的には $T \times X_t$ は $T$ の 1 倍（定数項 0 の線形結合）で表現できるため、完全な線形依存が生じます。

ただし介入後 $X_t = 1$ から、

$Z = T$

となるので介入前 $X_t = 0$ では線形依存は発生しません。

2.2 rank落ちによる問題

回帰分析では、デザイン行列が $X$ として表されますが、もしある説明変数が他の変数の完全な線形結合で表される場合、デザイン行列はフルランクではなくなります。
具体的には、Equation 1 の場合、説明変数として

定数項 $1$
時間変数 $T$
介入ダミー $X_t$
交互作用項 $T \times X_t$

が含まれますが、介入後では $T \times X_t$ = $T$ となるため、 $T$ と $T \times X_t$ の両方が同じ情報を持ちます。

その結果、行列のランクが本来の列数よりも低くなり、逆行列を求められなくなります。
このrank落ちは、回帰係数が一意に推定できない（識別できない）状態を招き、推定の信頼性に深刻な影響を与えます。

3. 再パラメータ化による解決策

上記の問題を解決するため、交互作用項を介入開始時点 $T_i$ で中心化する再パラメータ化を行います。具体的には、次の恒等式を利用します。

$T \times X_t = (T-T_i)X_t + T_iX_t$

これを Equation 1 に代入すると、

$Y_t = \beta_0 + \beta_1 T + \beta_2 X_t + \beta_3 [(T-T_i)X_t + T_iX_t$ ]

となり、整理すると、

$Y_t = \beta_0 + \beta_1 T + (\beta_2 + \beta_3 T_i) X_t + \beta_3 (T-T_i)X_t$

となります。ここで新たに、

$\beta'_2 = \beta_2 + \beta_3 T_i,\quad \beta'_3 = \beta_3$

と定義すると、モデルは以下の Equation 2 として書き換えられます。

$Y_t = \beta_0 + \beta_1 T + \beta'_2 X_t + \beta'_3 (T-T_i)X_t$

この形式では、介入直後（ $T = T_i$ ）に

$(T-T_i)X_t = 0$

となるため、即時効果が正しく $\beta'_2$ により評価されます。さらに、中心化された交互作用項 $(T-T_i)X_t$ は $T$ とは独立した変動を持つため、行列がフルランクに近い状態となります。これにより、推定が一意に行えるようになり、ランク落ちの問題が解消されます。

シミュレーションで確認します。


    library(car)

# シミュレーションパラメータ
T <- 0:150             # 0から150までの時間
T_i <- 75             # 介入開始時刻
X <- ifelse(T >= T_i, 1, 0)  # 介入ダミー：T < T_i は 0、T >= T_i は 1

# Equation 1 の交互作用項: T * X
interaction1 <- T * X

# Equation 2 の交互作用項: (T - T_i) * X
interaction2 <- (T - T_i) * X

# Equation 1 のデザイン行列 (Intercept は lm() が自動で追加)
df_eq1 <- data.frame(
  T = T,
  X = X,
  Interaction = interaction1
)

# Equation 2 のデザイン行列
df_eq2 <- data.frame(
  T = T,
  X = X,
  Interaction = interaction2
)

# ダミーの応答変数を生成（例：正規乱数）
set.seed(123)
df_eq1$Y <- rnorm(nrow(df_eq1))
df_eq2$Y <- rnorm(nrow(df_eq2))

# 線形モデルのフィッティング
model_eq1 <- lm(Y ~ T + X + Interaction, data = df_eq1)
model_eq2 <- lm(Y ~ T + X + Interaction, data = df_eq2)

# VIF の計算
vif_eq1 <- vif(model_eq1)
vif_eq2 <- vif(model_eq2)

cat("Equation 1 の VIF:\n")
print(vif_eq1)

cat("\nEquation 2 の VIF:\n")
print(vif_eq2)

以上からEquation2ではVIFが改善されているのが分かると思います。

以上、Interrupted Time-Series Analysisにおける多重共線性について考察しました。

2025-01-13

【Kaggle】Harmful Brain Activity 解法

こんにちは。私は情報系の学部4年生（院進学予定）です。研究室では主に数理統計の理論について取り組んでいます。そんな私が、Kaggleで開催されていた脳波（EEG）データを扱う「HMSコンペ」で奇跡的に入賞しました。

そしてkaggle Expartになりました。ヤッター

f:id:konchu2:20250308010746j:image

当時の取り組みを振り返りながら解法をまとめたいと思います。

はじめに

背景・自己紹介
- 学部：情報系（B4 → 院進学予定）
- 研究室：数理統計（理論中心）
Kaggleでの基本的な戦略
- フルスクラッチは難しいので、公開ノートブックやディスカッションを参考にする。
- 少しだけ改変できそうな部分は、AIに助言をもらいながら試行錯誤。
- いろいろとアンサンブルして Submit → あとは祈る。
- こうした地道な方針でも、意外と良い結果が得られることがあるのでおススメです。

この記事は、私が試行錯誤したプロセスを整理する目的で書いています。

コンペ概要

タスク

EEGデータ（脳波）を用いて、「発作（Seizure）」や「異常な脳活動」を分類するモデルを作る。
EEGデータをスペクトログラム（音の周波数のようなものを可視化した画像）に変換し、その画像を分類する。

ラベルの種類（6クラス）

Seizure（発作）
LPD（片側性周期性放電）
GPD（全般性周期性放電）
LRDA（片側性リズミカルデルタ活動）
GRDA（全般性リズミカルデルタ活動）
Other（その他）

注目ポイント：専門家の投票数

データには「専門家の投票数」という信頼性の指標が付与されており、多いほどデータの質が高いと考えられる。
“Other”に票が集中しているデータは曖昧なケースの可能性が高い。
投票数の扱い方を工夫することが、このコンペでの成績を左右すると感じました。

データの扱い方

高品質データと低品質データの分割

「投票数」を基準にデータを2つに分けました。

高品質データ：投票数が10票以上
低品質データ：投票数が1～9票

高品質データのみで先にモデルをトレーニングし、その後、低品質データをどう活用するかを工夫しました。この単純な戦略でも、モデルのパフォーマンスが大きく向上しました。

自己学習（擬似ラベリング）の活用

「低品質データ」をそのまま捨てるのではなく、モデルに疑似ラベルを予測させ、その結果を再度モデルに学習させる**自己学習（擬似ラベリング）**を試しました。大まかな流れは以下のとおりです。

初期トレーニング
- 高品質データのみでモデルを学習。
擬似ラベル付与
- 学習済みモデルで低品質データに対して予測ラベルを付与。
擬似ラベルの改良
- 低品質データの擬似ラベルと元のラベル情報を組み合わせて、新しいデータセットを作成。
- 高品質データの50%にも同様の処理を行い、データ全体のバランスを整える。
再トレーニング
- 改良したデータセットを使ってモデルを再学習。
ファインチューニング
- 最後に高品質データで微調整して仕上げる。

「本当に精度が上がるのか？」と半信半疑でしたが、結果はなかなか良好でした。まるでモデルを“育成”する感覚で、愛着が湧きました。

スペクトログラムの拡張

EEGデータはそのままだと扱いづらいため、スペクトログラムに変換するのが一般的です。しかし、提供された10分間のスペクトログラムだけでは物足りなかったため、複数の時間スケールでスペクトログラムを作り直しました。

10分間（全体像を把握）
50秒間（やや細部を確認）
10秒間（さらに細部を精密に確認）

これら3種類を組み合わせて512×512ピクセルの1枚の画像にまとめることで、「大まかな流れ」と「細かい特徴」の両方を捉えられるようにしました。結果、モデルの性能が大きく向上しました。

モデル選択と組み合わせ

EfficientNetB0

特徴：軽量で高速、かつ精度もそこそこ高い。
メリット：データが多くなくても、比較的良い性能を発揮。

ResNet

特徴：より深い構造で複雑なパターンを捉えられる。
メリット：EfficientNetB0だけでは拾えない詳細な特徴を補完。

なぜ両方を使ったのか？

EfficientNetB0 → 全体の大まかな特徴を素早く捉える。
ResNet → より深いレベルでのパターン抽出が可能。

この2つをアンサンブルすることで、速さと正確さの両方を両立しました。

振り返りと今後の展望

データの質を見極める重要性
- 専門家の投票数を基準にデータを分割しただけで、予想以上にモデルの性能が上がった。
- 「ラベルの信頼度」が高いデータを優先的に学習させるのは、精度向上のカギになると実感。
擬似ラベリング（自己学習）の可能性
- 低品質データも捨てずに活用することで、モデルを“自分で賢くさせる”アプローチが効果的だった。
- データをただ使うのではなく、“育てる”という発想が楽しい。
モデルの使い分け
- EfficientNetB0とResNetの組み合わせにより、速度と複雑なパターン抽出の両立ができた。
- 今後は時系列の情報をより活かせるモデルなどを試してみたい。

まとめ

以上が私が脳波（EEG）データを扱うHMSコンペで取り組んだ内容です。データの信頼度を見極めつつ、捨てずに活かす工夫をすることで、初めてのKaggleメダルを獲得することができました。機械学習にはまだまだ不慣れで、プログラミングも得意とはいえませんが、公開ノートブックやAIの助けを借りながらでも結果を出すチャンスは十分あると感じています。