import numpy as np
from pykalman import KalmanFilter

np.random.seed(42)
true = np.cumsum(np.random.randn(200)*0.2) + 10
obs = true + np.random.randn(200)*3

kf = KalmanFilter(
    transition_matrices=1, observation_matrices=1,
    transition_covariance=0.04, observation_covariance=9,
    initial_state_mean=obs[0])
est = kf.filter(obs)[0].flatten()

rmse = lambda x: np.sqrt(((x - true)**2).mean())
print(f"観測値のRMSE:   {rmse(obs):.4f}")
print(f"推定状態のRMSE: {rmse(est):.4f}")

簡単に説明します。

ライブラリのインポート

• numpy：数値計算を行うためのライブラリ。ここでは乱数生成やRMSEの計算に使用しています。
• pykalman.KalmanFilter：カルマンフィルタを実装したクラスです。

データの生成

• true = np.cumsum(np.random.randn(200)*0.2) + 10：真の潜在状態です。ランダムウォーク（各時点で小さくランダムに変動する）として生成しています。この値は現実には観測できないものとして扱います。
• obs = true + np.random.randn(200)*3：観測値です。真の状態に標準偏差3の大きな観測ノイズが加わっています。真の状態の変動（標準偏差0.2）に比べてノイズが非常に大きく、観測値だけを見ても真の状態がわかりにくい状況です。

カルマンフィルタの設定

• transition_matrices=1：状態遷移行列。1なので「次の状態は現在の状態と同じ（＋プロセスノイズ）」というランダムウォークモデルです。
• observation_matrices=1：観測行列。1なので「観測値は状態そのもの（＋観測ノイズ）」を意味します。
• transition_covariance=0.04：プロセスノイズの分散。真の状態の変動幅（0.2の2乗＝0.04）に対応します。
• observation_covariance=9：観測ノイズの分散。観測ノイズ（3の2乗＝9）に対応します。
• initial_state_mean=obs[0]：初期状態の推定値として最初の観測値を使用しています。

フィルタリングの実行

• kf.filter(obs)[0]：各時点までの観測値を使って状態を逐次推定します。戻り値の1つ目が推定状態の平均です。
• rmse：真の状態との二乗平均平方根誤差を計算するラムダ関数です。

(A) 観測データからトレンドと季節性を分解している → ×

時系列分解にはseasonal_decomposeやSTLを使います。カルマンフィルタは状態空間モデルの枠組みで潜在状態を推定する手法です。ただし、状態空間モデルを適切に設計すれば、トレンドや季節成分を状態として推定することも可能です。

(B) 観測ノイズを含むデータから、背後にある潜在状態を推定している → ○

正解です。カルマンフィルタは、ノイズに汚染された観測から、直接は観測できない真の状態を最適に推定します。出力のRMSEの改善（2.9648→0.7715）がその効果を示しています。

(C) 観測データの異常値を検出して除去している → ×

異常値の検出には別の手法を使います。カルマンフィルタはすべての観測を使って状態を推定するもので、異常値の除去を目的としたものではありません。

(D) 観測データの周波数成分を分析している → ×

周波数成分の分析にはFFTやWelch法などのスペクトル解析を使います。カルマンフィルタは時間領域で逐次的に状態を推定する手法です。

カルマンフィルタは、状態空間モデルという枠組みの上で動作します。状態空間モデルは2つの方程式から構成されます。

状態方程式（システムモデル）:
$$
x_t = F x_{t-1} + w_t, \quad w_t \sim N(0, Q)
$$

潜在状態x_tが時間とともにどう変化するかを記述します。Fは状態遷移行列（コードのtransition_matrices）、Qはプロセスノイズの共分散（transition_covariance）です。

観測方程式（観測モデル）:

$$
y_t = H x_t + v_t, \quad v_t \sim N(0, R)
$$

潜在状態x_tがどのように観測値y_tとして現れるかを記述します。Hは観測行列（observation_matrices）、Rは観測ノイズの共分散（observation_covariance）です。

今回のコードではF = 1、H = 1としているため、「状態はランダムウォークし、観測はその状態にノイズが乗ったもの」という最もシンプルなモデルになっています。

カルマンフィルタは、各時点で「予測」と「更新」の2段階を繰り返します。

1. 予測ステップ

前時点までの情報から、現時点の状態を予測します。

$$
\hat{x}_{t|t-1} = F \hat{x}_{t-1|t-1}
$$

同時に、予測の不確実性（誤差共分散）も更新されます。時間が経過するほど不確実性は増大します。

2. 更新ステップ

新しい観測値y_tが得られたら、予測を修正します。

$$
\hat{x}_{t|t} = \hat{x}_{t|t-1} + K_t (y_t – H \hat{x}_{t|t-1})
$$

ここでK_tはカルマンゲインと呼ばれ、「予測と観測のどちらをどれだけ信頼するか」を決める重みです。

• 観測ノイズRが大きい（観測が不正確） → カルマンゲインは小さくなり、予測を重視する。
• プロセスノイズQが大きい（状態が不安定） → カルマンゲインは大きくなり、観測を重視する。

今回のコードではR = 9とQ = 0.04で、観測ノイズが圧倒的に大きい設定です。そのためカルマンゲインは小さくなり、観測値の変動に振り回されず滑らかな推定が得られています。

ノイズ除去・平滑化

移動平均や指数平滑法と同様に平滑化に使えますが、カルマンフィルタはノイズの統計的性質をモデルに組み込むため、より理論的に最適な推定が得られます。

欠損値の補間

観測が欠損している時点でも、予測ステップのみを実行することで状態の推定を継続できます。欠損値に頑健であることは、カルマンフィルタの実用上の大きな利点です。

時変パラメータの推定

回帰係数などのパラメータが時間とともに変化する場合、そのパラメータを状態として扱うことで、時変係数を逐次推定できます。

ARIMAとの関係

statsmodelsのARIMAやSARIMAXは、内部的に状態空間表現に変換され、カルマンフィルタで尤度が計算されています。つまり、これまでのノックで扱ってきたARIMAモデルの推定も、裏側ではカルマンフィルタが動いているのです。

フィルタリングとスムージング

kf.filter()は各時点までの情報のみを使う「フィルタリング」で、リアルタイム処理に適しています。一方、kf.smooth()は全期間のデータを使って過去の状態を推定し直す「スムージング」で、事後分析ではこちらの方が高精度な推定が得られます。