40 人が閲覧(直近 30 日) Python

カルマンフィルタとは?仕組みを1次元の例でわかりやすく解説【Python実装とQ・Rの決め方】

カルマンフィルタとは?仕組みを1次元の例でわかりやすく解説【Python実装とQ・Rの決め方】

カルマンフィルタは、「モデルで予測した値」と「センサーで測った値」を、それぞれの不確かさ(分散)に応じた比率で混ぜて、見えない状態を1ステップずつ推定する手法です。GPSの位置飛びをならす、速度を測っていないのに位置の列から速度を推定する、といった場面で使われます。

この記事では、まず温度計1本の数値例で「なぜその比率で混ぜるのか」を示し、次に行列の式、NumPyで30行ほどの実装、実測したQとRの効き方の順に進めます。位置推定の単一試行は記事中のコードで再現できます。移動平均比較、Q・Rの変更、急加速、外れ値の各実験は、別の検証スクリプトで計算した結果です。

まとめ:カルマンフィルタの要点

  • カルマンフィルタは「予測」と「更新」の2ステップを繰り返す逐次推定法です。予測ではモデルに従って共分散を伝播させてQを加え、更新では観測情報を反映して共分散を減らします。
  • 混ぜる比率がカルマンゲインKです。予測の分散が4、観測の分散が1ならK=0.8で、推定値は観測側へ8割寄ります。
  • 線形モデルとガウス雑音の条件では平均二乗誤差が最小の推定になり、雑音がガウスでなくても線形推定の中では最良です。
  • 実測では、位置を標準偏差5mで観測する等速度モデルで位置RMSEが4.99mから2.73mに下がりました。窓5点の移動平均は6.41mと、観測そのままより悪化しました。
  • QとRの相対的な大きさはゲインを左右しますが、初期過渡では初期共分散P0も影響します。絶対値が妥当かどうかは、イノベーションを正規化した値(NIS)の平均が観測の次元数に近いかで点検します。
  • 記録済みのデータを後からまとめて処理するなら、未来の観測も使うRTS平滑化が使えます。同じ設定で位置RMSEは2.73mから1.50mに下がりました。QとRが分からない時系列は、statsmodelsで最尤推定できます。
  • 非線形なら拡張カルマンフィルタ(EKF)かUKF、分布が多峰なら粒子フィルタを選びます。

カルマンフィルタとは:予測と観測を不確かさで重み付けする逐次推定

カルマンフィルタが扱うのは、「状態は直接見えないが、時間とともにどう変わるかのモデル」と「ノイズの乗った観測」がある問題です。毎時刻、モデルで次の状態を予測し、観測が来たら予測とのずれを使って修正します。過去の観測をすべて保存しておく必要はなく、手元に持つのは現在の推定値と、その不確かさを表す共分散行列だけです。

最適性の条件も明確です。NASA NTRSに収録されたKimらのIEEE Access論文は、ガウス白色雑音を受ける離散時間の線形システムでカルマンフィルタが平均二乗誤差を最小にする推定器(MMSE推定器)であり、雑音がガウスでなくても平均と共分散が既知なら「線形推定アルゴリズムの中で最良」だと整理しています(NASA NTRS 20230017116)。確率の見方では、予測を事前分布、観測を尤度として事後分布を求めるベイズの定理を、正規分布に限って閉じた式で解いたものに当たります。

最小二乗法との違いもここから説明できます。一定のスカラー量を直接観測するF=H=1、Q=0のモデルで、観測誤差が無相関かつ等分散であり、事前情報を持たない極限では、推定値は標本平均、つまり最小二乗解と一致します。事前の分散を10⁹として50点の観測で試すと、推定値と標本平均の差は各時点で最大2.3×10⁻⁸、ゲインと1/nの差は最大4×10⁻⁹でした。事前の分散が有限なので、厳密には近似的な一致です。カルマンフィルタは最小二乗法を「状態が動く」「観測が逐次届く」場合へ広げたものと捉えると分かりやすくなります。

1960年のKalman論文とApollo計画での実用化

名前の由来はルドルフ・カルマン(Rudolf E. Kálmán、1930〜2016年)です。原論文は「A New Approach to Linear Filtering and Prediction Problems」で、Journal of Basic Engineering 第82巻第1号(1960年3月)35〜45頁に掲載されました。カルマンはこの業績などで2008年の米国国家科学賞(National Medal of Science)と同年のドレイパー賞を受けています。

最初に公に知られた応用は、1960年代初頭のNASAエイムズ研究センターによるApollo宇宙船の月周回航法の研究です。NASAの技術メモ(McGee・Schmidt、NASA-TM-86847、1985年)によると、Stanley F. Schmidtらは非線形な軌道の式を、推定中の現在の状態のまわりで線形化し直しながら使う方式を作り、これが後に「拡張カルマンフィルタ」と呼ばれるようになりました。つまりEKFは理論の後付けではなく、実用化の最初から一緒に生まれています。

カルマンフィルタの主な応用先と推定する状態

使いどころは「観測できる量」と「知りたい量」が違う場面です。位置推定の例では、GPSで位置を、IMUで加速度と角速度を測り、両者を統合して位置・速度・姿勢を推定します。ロボットが地図を作りながら自己位置を推定するSLAMでも、初期の代表的な方式はEKFを使うものでした。

分野 観測する量 推定する状態
車両・ドローンの航法 GPS位置、IMUの加速度・角速度 位置・速度・姿勢
ロボットの自己位置推定 車輪の回転数、LiDARやカメラの特徴点 自己位置と地図の目印の位置
バッテリー管理 端子電圧・電流 充電状態(SOC)
画像の物体追跡 フレームごとの検出位置 位置・速度(見失った間の予測)
経済・需要の時系列 月次の実績値 トレンド・季節成分

システムとして導入する場合は、カルマンフィルタの実装に加え、センサーからのデータ収集・時刻合わせ・推定結果の配信も開発範囲に含めて見積もります。こうした収集から推定・可視化までの開発は、弊社のAI/IoTソリューションでも支援しています。

1次元の温度推定:予測値と観測値の重み付け

部屋の温度を推定する場面で考えます。前回の推定とモデルから、今の温度は20.0℃と予測しました。ただしこの予測には標準偏差2℃ぶんの自信のなさがあり、分散は4です。そこへ温度計が21.0℃を示しました。温度計の誤差は標準偏差1℃、分散は1とします。

どちらを信じるかを、分散の比で決めます。カルマンゲインは K = 予測の分散 / (予測の分散 + 観測の分散) = 4 / (4 + 1) = 0.8 です。推定値は予測から観測へ向けてずれの8割だけ動かし、20.0 + 0.8 × (21.0 − 20.0) = 20.8℃ になります。混ぜた後の分散は (1 − 0.8) × 4 = 0.8 で、予測の4とも観測の1とも比べて小さくなります。2つの不確かな情報を合わせると、どちらか一方より確かな値が得られる点が要です。観測の分散が予測より大きい場合はKが0.5未満になり、推定値は予測寄りに留まります。

重み0.8の導出:推定誤差の分散最小化

0.8という比率は、混ぜた結果の分散が最も小さくなる重みとして導けます。予測に重み (1 − w)、観測に重み w をかけて足すとします。予測の誤差と観測の誤差が無相関なら、混ぜた値の誤差の分散は (1 − w)² × 4 + w² × 1 です。これをwで微分して0とおくと −8(1 − w) + 2w = 0 から w = 4 / (4 + 1) = 0.8 となり、カルマンゲインと一致します。そのときの分散は 0.2² × 4 + 0.8² × 1 = 0.8 で、前段の値と同じです。

比較のため、重みを0.5にして単純に平均すると、分散は 0.25 × 4 + 0.25 × 1 = 1.25 になり、温度計1本の分散1より悪くなります。不確かさを無視して混ぜると、良い方の情報まで薄めてしまうということです。予測の分散をP、観測の分散をRと書けば最適な重みは P / (P + R) で、行列版のカルマンゲインはこれを多次元に広げたものです。

4ステップの数値の流れ:ゲインが下がっていく様子

次の時刻には、この20.8℃(分散0.8)から再び予測します。時間が経つ間に温度が変わる可能性があるので、予測の分散には「モデルの不確かさ」Q(ここでは0.1)を足します。温度は一定とみなして予測値は前回の推定値のまま、観測値を21.0、20.6、21.3、20.9℃とした4ステップの計算が次の表です。

ステップ 予測値 予測の分散 ゲインK 観測値 推定値 推定の分散
1 20.00 4.000 0.800 21.0 20.80 0.800
2 20.80 0.900 0.474 20.6 20.71 0.474
3 20.71 0.574 0.365 21.3 20.92 0.365
4 20.92 0.465 0.317 20.9 20.92 0.317

予測で分散が増え(0.800→0.900)、観測で縮む(0.900→0.474)。この繰り返しがカルマンフィルタの全体です。推定の分散が小さくなるにつれてゲインも0.800、0.474、0.365、0.317と下がり、観測1本が推定を動かす力は弱まっていきます。4ステップ目は観測20.9℃が予測20.92℃とほぼ同じなので、推定値は表示桁では動きません。なお、この例はR=1なので推定の分散がKと同じ数値になります((1 − K) × 予測の分散 = K × R が成り立つため)。

予測ステップと更新ステップの式:状態空間モデルの行列版

1次元の例を、位置と速度のような複数の量へ広げたものが行列版です。記号が増えるだけで、やっていることは前章と同じです。

状態方程式と観測方程式:F・H・Q・Rの意味

カルマンフィルタは、系を2本の式で表す状態空間モデルを前提にします。状態方程式 x_k = F x_(k−1) + w_k は「状態が1ステップでどう変わるか」、観測方程式 z_k = H x_k + v_k は「状態のどの部分がどう見えるか」を表します。w_kとv_kは平均0の白色雑音で、それぞれの共分散がQとRです。ここで示す式では、両雑音と初期状態の誤差は互いに無相関とします。

記号 名前 等速度モデルでの中身
x 状態ベクトル [位置, 速度]
F 状態遷移行列 位置 += 速度 × dt
H 観測行列 [1, 0](位置だけ測る)
Q プロセスノイズ共分散 加速度のゆらぎの大きさ
R 観測ノイズ共分散 センサー誤差の分散(5mの2乗=25)
P 誤差共分散行列 推定値の不確かさ(毎回更新)

観測行列Hが「速度を測っていない」ことを表している点に注目してください。それでも推定値には速度が入ります。位置の観測が続くと、位置の変化から速度が間接的に絞り込まれるためで、これは状態方程式で位置と速度をつないでいるから可能になります。

5本の更新式とカルマンゲインの意味

1ステップの計算は次の5本です。添え字 pred は予測後、無印は更新後の値を表します。

段階 式 意味
予測 x_pred = F x 状態を1ステップ進める
予測 P_pred = F P Fᵀ + Q 不確かさをFで運びQを足す
更新 K = P_pred Hᵀ (H P_pred Hᵀ + R)⁻¹ 混ぜる比率を決める
更新 x = x_pred + K (z − H x_pred) 予測を観測側へ寄せる
更新 P = (I − K H) P_pred 不確かさが縮む

3本目のカルマンゲインは、1次元の 予測の分散 / (予測の分散 + 観測の分散) を行列に置き換えたものです。括弧の中の S = H P_pred Hᵀ + R はイノベーション(観測と予測の差 z − H x_pred)の分散で、後の章でQとRの点検に使います。H=1の1次元モデルでは、他の条件を固定すると、予測分散が大きいほどKが大きくなり、観測分散Rが大きいほどKが小さくなります。多次元では、観測行列と共分散の成分間の関係もゲインに影響します。

PythonとNumPyで実装する等速度モデルの位置推定

直線上を動く対象の位置を、標準偏差5mの誤差で毎ステップ観測する例を実装します。速度は一定ではなく、ランダムな加速度でゆっくり変わる設定です。カルマンフィルタ専用ライブラリは使わず、NumPyで前章の5本の式を書きます。

import numpy as np

dt, q, r = 1.0, 0.1, 25.0                     # 時間刻み・加速度ノイズの強さ・観測ノイズ分散(5mの2乗)
F = np.array([[1.0, dt], [0.0, 1.0]])         # 状態遷移行列:位置 += 速度*dt
H = np.array([[1.0, 0.0]])                    # 観測行列:位置だけを測る
Q = q * np.array([[dt**3 / 3, dt**2 / 2], [dt**2 / 2, dt]])
R = np.array([[r]])

def kalman_filter(zs, x, P):
    I = np.eye(len(x))
    est, gains = [], []
    for z in zs:
        x = F @ x                             # 予測:状態を1ステップ進める
        P = F @ P @ F.T + Q                   # 予測:不確かさを F で運び Q を足す
        S = H @ P @ H.T + R                   # イノベーション(観測-予測)の分散
        K = P @ H.T @ np.linalg.inv(S)        # カルマンゲイン
        x = x + K @ (z - H @ x)               # 更新:予測を観測側へ K だけ寄せる
        A = I - K @ H
        P = A @ P @ A.T + K @ R @ K.T         # 更新:不確かさが縮む(Joseph形式)
        est.append(x.copy()); gains.append(K[0, 0])
    return np.array(est), np.array(gains)

# 真値と観測を作る(等速度+ランダムな加速度、位置を標準偏差5mで観測)
rng = np.random.default_rng(0)
x_true, xs, zs = np.array([0.0, 1.0]), [], []
for _ in range(200):
    x_true = F @ x_true + rng.multivariate_normal([0, 0], Q)
    xs.append(x_true.copy())
    zs.append(H @ x_true + rng.normal(0, np.sqrt(r), 1))
xs, zs = np.array(xs), np.array(zs)

est, gains = kalman_filter(zs, x=np.zeros(2), P=np.diag([100.0, 10.0]))
rmse = lambda a, b: np.sqrt(np.mean((a - b) ** 2))
print(f"観測そのまま  位置RMSE {rmse(zs[:, 0], xs[:, 0]):.2f} m")
print(f"カルマン推定  位置RMSE {rmse(est[:, 0], xs[:, 0]):.2f} m  速度RMSE {rmse(est[:, 1], xs[:, 1]):.2f} m/step")
print("ゲインの推移", np.round(gains[[0, 1, 4, 9, 19, 199]], 3))

Python 3.9.6とNumPy 2.0.2で実行した結果です。

観測そのまま  位置RMSE 4.88 m
カルマン推定  位置RMSE 2.53 m  速度RMSE 0.65 m/step
ゲインの推移 [0.815 0.572 0.509 0.354 0.3   0.299]

位置の誤差はほぼ半分になり、観測していない速度も誤差0.65m/stepで推定できています。共分散の更新は、表の5本目 P = (I − K H) P_pred と数学的に同じで、丸め誤差でPの対称性が崩れにくいJoseph形式で書いています。同じ入力をfilterpy 1.4.5の KalmanFilter に与えると、全200ステップで推定値の差は0でした。

ゲインが一定値に収まる:定常カルマンフィルタ

出力の最終行を見ると、ゲインは最初の0.815から20ステップ目に0.300まで下がり、その後は表示した小数第3位まででは0.299に落ち着きます。初期値の不確かさ(P0の対角に100と10)が観測で消え、FやQ、Rが時間で変わらない系では、予測で増える分散と観測で減る分散が釣り合うためです。

F・H・Q・Rが一定で、可検出性など定常解への収束条件を満たす場合は、収束後のゲインを事前に計算しておき、行列の逆行列を毎回解かずに固定の係数で更新する実装が取れます。マイコンのように計算資源が限られる環境でよく使われる形で、位置と速度に固定ゲインをかけるα-βフィルタは、この定常解に相当します。

移動平均との比較:移動対象の追従遅れと位置誤差

ノイズを減らすだけなら移動平均でもよいのでは、という疑問に数値で答えます。上と同じ設定で乱数の種を500通り変え、最初の20ステップを除いた位置RMSEを平均しました。移動平均は時刻kまでの直近w点の平均(後ろ向きの窓)です。再現するには、上のコードの観測生成を default_rng(seed) の種0〜499で回し、各回で kalman_filter と移動平均のRMSEを21ステップ目以降について計算して平均します。

手法 位置RMSE(500試行平均)
観測そのまま 4.99 m
カルマンフィルタ 2.73 m
移動平均(窓5点) 6.41 m
移動平均(窓10点) 13.26 m
移動平均(窓20点) 27.10 m

移動平均は窓を広げるほどノイズは消えますが、動いている対象に対しては窓幅のおよそ半分の時間ぶん古い位置を返すため、その遅れがそのまま誤差になります。窓5点ですでに観測そのままより悪化しました。カルマンフィルタは速度も状態として持ち、「今はこのあたりまで進んでいるはず」と予測してから混ぜるので、この実験では後ろ向き移動平均より小さい誤差で現在位置を推定できました。ただし、モデルが外れると追従遅れも生じます。移動平均との選択は、許容する遅れとモデル化の可否に応じて判断します。

記録済みデータのRTS平滑化:後続観測による推定誤差の低減

カルマンフィルタは時刻kの推定に、時刻kまでの観測しか使いません。リアルタイム処理ではそれが前提ですが、走行ログや過去の計測値を後から解析する場合は、時刻kより後の観測も使えます。そのための代表的な手法が、Rauch・Tung・Striebelが1965年に発表したRTS平滑化(固定区間平滑化)です(AIAA Journal 第3巻第8号)。まず通常どおり最後までフィルタを回し、各時刻の予測値・推定値と共分散を保存します。次に最終時刻から逆向きに、C = P_k Fᵀ P_pred(k+1)⁻¹ を係数として、1つ後の時刻の平滑化結果との差で推定値を補正していきます。

前述の等速度モデルで、乱数の種500通り・21ステップ目以降について測った結果です。

手法 位置RMSE 速度RMSE
カルマンフィルタ 2.73 m 0.72 m/step
RTS平滑化 1.50 m 0.38 m/step

位置の誤差は約45%、速度の誤差は約47%減りました。種0の試行で101ステップ目(配列添え字100)の位置の分散を見ると、フィルタの7.48に対して平滑化は2.22です。一方、最終時刻ではその後の観測が無いので、平滑化の結果はフィルタと一致しました。平滑化の改善は系列の中ほどで大きく、末端では小さくなります。リアルタイムの制御には使えないので、オンラインはフィルタ、オフラインの解析は平滑化と使い分けます。pykalmanの smooth() やstatsmodelsの smoothed_state で、同じ処理を自分で書かずに使えます。

QとRの設定:相対的な大きさとNISによる整合性の点検

カルマンフィルタの実装で最も迷うのがQとRの値です。Rは、静止状態で取った観測の分散などから見積もります。仕様書を使う場合は、精度の表記が標準偏差・信頼区間・最大誤差のどれかを確認し、ランダムな観測誤差の分散へ換算します。一定の偏りは別に補正またはモデル化します。難しいのは、モデルがどれだけ外れるかを表すQの方です。そこで、真のQ(q=0.1)とRを知っている上の設定で、フィルタに渡すQだけを何倍かにずらして500試行ずつ測りました。

フィルタに渡したQ 定常ゲイン 位置RMSE NIS平均
真値の0.01倍 0.106 7.47 m 3.58
真値の0.1倍 0.181 3.63 m 1.44
真値(1倍) 0.299 2.73 m 1.00
真値の10倍 0.469 3.13 m 0.86
真値の100倍 0.675 3.81 m 0.72

Qを真値に合わせたときにRMSEが最小になり、小さすぎる側の悪化の方が急です。もう一つ分かったのは、Qを10倍にした場合とRを0.1倍にした場合で、定常ゲイン0.469・位置RMSE3.13mが表示桁で一致したことです。初期状態とモデル・観測列を固定し、Q・R・初期共分散P0をすべて同じ正の倍率で変えれば、推定値の軌跡は変わりません。P0を固定したままQとRだけを変えると、初期過渡の推定値は変わります。上の比較でも、同じ観測列での推定値の差は最初の20ステップで最大1.28m、21ステップ目以降は0.0004m以下でした。

ただし、フィルタが自分で報告する不確かさPは絶対値に比例して変わります。これを点検するのがNIS(正規化イノベーション二乗)で、毎ステップ yᵀ S⁻¹ y(yはイノベーション)を計算して平均します。QとRが正しければ平均は観測の次元数(ここでは1)に近づき、実測でも1.00でした。同じ比率のQ10倍とR0.1倍でも、NISは0.86と8.55で大きく違います。NIS平均が1から継続的に外れる場合は、標本数によるばらつきに加え、モデル誤差・バイアス・外れ値も点検します。QやRの変更方向は、これらの原因と他の検証指標を確認して決めます。

実務では、Rを実測分散で見積もり、Qを変えてNISを点検する方法が出発点になります。ただし平均の一致だけでは調整の正しさを保証できず、イノベーションの平均・自己相関や、別の検証データでの推定誤差も確認します。

QとRをデータから推定する:statsmodelsの最尤推定

観測の系列があれば、QとRを手で探さずにデータから推定する方法もあります。カルマンフィルタは毎ステップのイノベーションとその分散Sを出すので、観測列全体の尤度(そのQ・Rのもとで観測列が得られる確からしさ)を計算できます。尤度が最大になるQ・Rを探すのが最尤推定です。statsmodelsの公式ドキュメントも、UnobservedComponents で追加で最尤推定するのは確率的な成分の分散だと説明しています(statsmodels: UnobservedComponents)。

水準がランダムウォークで動き、そこに雑音が乗る「ローカルレベルモデル」で試します。真の値は、水準のゆらぎの分散(Qに当たる)を0.5、観測雑音の分散(R)を4.0としました。

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(1)
level = np.cumsum(rng.normal(0, np.sqrt(0.5), 500))   # 真の水準:分散0.5で揺れるランダムウォーク
y = level + rng.normal(0, 2.0, 500)                     # 観測:分散4.0の雑音を足す

res = sm.tsa.UnobservedComponents(y, level="local level").fit(disp=False)
print({k: round(float(v), 3) for k, v in zip(res.model.param_names, res.params)})
print("観測そのまま", np.sqrt(np.mean((y - level) ** 2)).round(2),
      "フィルタ", np.sqrt(np.mean((res.filtered_state[0] - level) ** 2)).round(2),
      "平滑化", np.sqrt(np.mean((res.smoothed_state[0] - level) ** 2)).round(2))

statsmodels 0.15.0(Python 3.12.13)と0.14.6(Python 3.9.6)のどちらでも、出力は同じでした。

{'sigma2.irregular': 4.431, 'sigma2.level': 0.33}
観測そのまま 2.11 フィルタ 1.04 平滑化 0.83

sigma2.irregular がR、sigma2.level がQの推定値です。観測500点の1系列では、Rは真値4.0に対して4.431、Qは0.5に対して0.33でした。乱数の種を50通り変えて推定し直すと、Rの推定値は3.348〜4.604(平均3.992)、Qは0.233〜0.734(平均0.489)に散らばりました。平均はどちらも真値に近い一方、Qは1系列だけだと真値の半分から1.5倍程度まで振れます。最尤推定の値をそのまま固定せず、上のNISや、次項のように急な動きを含むデータで確かめてから使うのが安全です。時系列分析の手法全体の中での位置づけは時系列分析の手法の使い分けで扱っています。

Qの過小設定による失敗:急な速度変化への追従遅れ

Qを小さくすると平時の推定は滑らかになり、一見よく見えます。そこで、速度が1m/stepで100ステップ進んだ後、突然5m/stepに上がる動きを200試行で与えました。

フィルタに渡したQ 平時のRMSE(20〜99) 加速後30ステップの最大絶対誤差(200試行平均)
基準の0.01倍 1.45 m 22.17 m
基準(1倍) 2.37 m 7.85 m
基準の100倍 3.77 m 8.68 m

Qを基準の0.01倍にしたフィルタは、平時の誤差が最も小さい一方で、加速後30ステップの最大絶対誤差は200試行平均で22.17mでした。モデルを信じすぎて、観測が「速くなった」と告げても予測側に留まるためです。平時のデータだけでQを詰めると、この小さすぎる側に寄りがちです。チューニングには、想定する最も急な動きを含むデータを必ず入れてください。

非線形のときの選択肢:拡張カルマンフィルタ・UKF・粒子フィルタ

ここまでの式は、状態の変化と観測がどちらも行列の掛け算で書ける線形の場合です。ロボットの向きを含む運動や、距離と角度で測るレーダーのように、状態や観測に三角関数が入ると線形ではなくなります。その場合の選択肢は3つです。

手法 非線形の扱い方 向く場面
拡張カルマンフィルタ(EKF) 現在の推定値のまわりでヤコビ行列により線形化 非線形が弱く、微分が書ける系
UKF(無香料カルマンフィルタ) 少数のシグマ点を非線形関数に通して平均・共分散を求める 非線形が強い、微分が面倒な系
粒子フィルタ 多数のサンプル(粒子)で分布そのものを表す 分布が多峰・非ガウス

EKFは前述のとおりApolloの航法で生まれた方式で、共分散の計算ではFとHをヤコビ行列に置き換えますが、状態の予測には非線形関数fを、イノベーションには非線形観測関数hを使います。弱点は線形化の誤差で、初期値が真値から遠いと発散することがあります。UKFはJulierとUhlmannが1997年に発表した方式で、微分を使わず、代表的な構成ではn次元に対して2n+1個のシグマ点を使います。ただし、雑音を状態に含める拡大状態方式では、その拡大後の次元に応じて点数が増えます。「無香料カルマンフィルタ」という呼び方は、Unscented Kalman Filterの直訳です。

粒子フィルタはモンテカルロ法で分布を近似する手法で、「対象は通路Aか通路Bのどちらかにいる」のような二つ山の分布を扱えます。代わりに粒子数に比例して計算が増えます。どれを選ぶかは、まず線形で書けるかを確かめ、書けなければEKF、EKFで発散するならUKF、分布が一つの山で表せないときに粒子フィルタ、の順で検討するのが無駄のない進め方です。

ライブラリの選び方:filterpy・pykalman・OpenCV・statsmodels・MATLAB

前章までの実装で分かるとおり、線形カルマンフィルタ自体は数十行で書けます。ライブラリを選ぶ基準は「フィルタ以外に何を任せたいか」です。版と日付は2026年10月時点の各公式ページ(PyPI・GitHubのリリース)で確認した値です。

ライブラリ 最新版(公開日) 任せられること
filterpy 1.4.5(2018-10-10) KF・EKF・UKFの学習用実装
pykalman 0.11.2(2026-01-31) 平滑化、EMによるパラメータ推定
OpenCV 5.0.0(2026-06-06) 画像処理と組み合わせた物体追跡
statsmodels 0.15.0(2026-08-27) 時系列モデルのノイズ分散を最尤推定
MATLAB R2026b 制御系設計からの定常ゲイン計算

filterpyは解説書と対応していて読みやすい一方、PyPIの最新版は2018年10月の1.4.5のままです。新規のプロダクトで長く保守する場合は、リリース日だけでなく、依存関係への対応、テスト、修正を担う体制まで比較して選んでください。OpenCVの cv2.KalmanFilter(dynamParams, measureParams, controlParams, type) は、5.0.0でも4系最新の4.14.0(2026-07-19)でも同じ引数で、predict() と correct(measurement) の2メソッドで上記の5本の式を回します。既定の型は32bit浮動小数点(CV_32F)なので、値の桁が大きい用途では型を指定してください。

QとRが分からないという悩みに直接答えるのはstatsmodelsで、前述のとおりノイズの分散の最尤推定と平滑化までを1つのモデルで扱えます。MATLABでは、Control System Toolboxの kalman 関数がプラントモデルとQ・Rからフィルタを設計します。制御系の開発環境としての全体像はMATLAB/Simulinkの解説にまとめています。

カルマンフィルタが向かない場面と外れ値への備え

カルマンフィルタは前提が崩れると、何も警告を出さずに間違った値を返します。採用を見送るか、手当てが要るのは次の場面です。

  • 外れ値が混ざる観測:通常の更新式は観測と予測の差にゲインを掛けて補正するため、想定した雑音の大きさを超える外れ値にも推定値が引っ張られます。後述のゲートで対処します。
  • 状態の変化をモデルにできない:状態変化を妥当なモデルで表せない場合、推定精度を保証できません。ゲインは設定したQ・Rから計算されるだけなので、モデルが外れていてもフィルタは気づかず、外れた予測に引っ張られた値を返し続けます。
  • 分布が二つ以上の山を持つ:位置の候補が複数ある状況は正規分布1つで表せません。粒子フィルタを使います。
  • 観測の時刻がずれる:センサーごとに遅延が違うのに受信順で更新すると、古い観測を新しい予測に混ぜることになります。タイムスタンプで並べ替えるか、遅延を扱う拡張が必要です。

外れ値の影響を測りました。前述の設定で151ステップ目(配列添え字150)の観測だけに+100mの外れ値を足すと、その時点の推定誤差は500試行平均で29.8mになりました。定常ゲインが約0.3なので、外れ値の3割がそのまま推定に入った計算です。

対策は、NISを使った観測のゲート(棄却)です。更新の前にイノベーションの yᵀ S⁻¹ y を計算し、自由度1のカイ二乗分布の99.9%点10.83を超えたら、その観測は捨てて予測値を使います。この処理を入れると同じ時点の誤差は2.6mまで下がりました。200ステップ中に棄却された観測は平均1.2個で、外れ値1個と、正常な観測の誤棄却(0.1%の確率)がおおむね0.2個です。センサーデータの異常を別の仕組みで検出する方法は時系列データの異常検知で扱っています。

カルマンフィルタに関するよくある質問

カルマンフィルタと最小二乗法の違いは何ですか?

最小二乗法は残差の二乗和を最小にする推定基準で、逐次計算する方式もあります。カルマンフィルタは状態の時間変化と雑音をモデル化し、観測ごとに推定を更新します。F=H=1、Q=0で一定のスカラー量を直接観測し、観測誤差が無相関かつ等分散なら、事前情報を持たない極限で最小二乗解(標本平均)と一致します。状態が動く場合に広げたものがカルマンフィルタです。

状態空間モデルとカルマンフィルタの関係は何ですか?

状態空間モデルは、系を状態方程式 x_k = F x_(k−1) + w_k と観測方程式 z_k = H x_k + v_k の2本で書く「モデル」です。カルマンフィルタは、そのモデルが線形で雑音がガウスのときに、観測から状態を逐次推定する「計算手順」に当たります。同じ計算で観測列の尤度も求まるため、状態空間モデルのノイズの分散を最尤推定するときにも使われます。

カルマンフィルタと拡張カルマンフィルタの違いは何ですか?

カルマンフィルタは状態の変化と観測が線形の場合の手法で、拡張カルマンフィルタ(EKF)は非線形の式を現在の推定値のまわりで線形化して同じ計算を使います。EKFはApollo計画の航法研究の中で生まれました。

カルマンゲインとは何ですか?

観測と予測の差を状態の修正量へ変換する係数です。状態を同じ単位で直接観測する1次元モデル(H=1)では 予測の分散 / (予測の分散 + 観測の分散) で、0〜1の値を取ります。1に近いほど観測を、0に近いほど予測を信じます。

Pythonでカルマンフィルタを使うにはどうすればよいですか?

線形ならNumPyで数十行で書けます。ライブラリならfilterpy(最新版は2018年の1.4.5)、pykalman、OpenCVの cv2.KalmanFilter、ノイズ分散も推定したい時系列ならstatsmodelsが選択肢です。

関連記事

お気に入りに入れた記事の一覧

この記事は以下の記事からリンクされています

資料請求

今日のトレンド記事 直近 24 時間で、いつもより多く読まれている記事

  1. 2026.10.23 テックブログ Renta!不正アクセスでサービス停止:再開未定の現状と運営側の停止判断・証拠保全・段階的再開
  2. 2026.07.09 テックブログ Muse Imageとは?Metaの画像生成AIの使い方・API料金(1枚0.01ドル)と日本での提供状況【2026年10月】
  3. 2026.10.22 テックブログ Tomcatのバージョン一覧とサポート期限【2026年10月】最新版11.0.27、9.0系は2027年3月末終了・9.1系は2030年末まで
  4. 2025.12.10 コラム 今の50代は何世代?バブル世代・団塊ジュニア・氷河期の境目と特徴
  5. 2026.10.24 テックブログ ブラストメール不正アクセス:配信先288件の漏えいと利用企業の報告・通知の手順

RELATED POSTS 関連記事

目次