カルマンフィルタは、「モデルで予測した値」と「センサーで測った値」を、それぞれの不確かさ(分散)に応じた比率で混ぜて、見えない状態を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)の平均が観測の次元数に近いかで点検します。
- 非線形なら拡張カルマンフィルタ(EKF)かUKF、分布が多峰なら粒子フィルタを選びます。
カルマンフィルタとは:予測と観測を不確かさで重み付けする逐次推定
カルマンフィルタが扱うのは、「状態は直接見えないが、時間とともにどう変わるかのモデル」と「ノイズの乗った観測」がある問題です。毎時刻、モデルで次の状態を予測し、観測が来たら予測とのずれを使って修正します。過去の観測をすべて保存しておく必要はなく、手元に持つのは現在の推定値と、その不確かさを表す共分散行列だけです。
最適性の条件も明確です。NASAの技術文書は、ガウス白色雑音を受ける離散時間の線形システムでカルマンフィルタが平均二乗誤差を最小にする推定器(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) |
| 画像の物体追跡 | フレームごとの検出位置 | 位置・速度(見失った間の予測) |
| 経済・需要の時系列 | 月次の実績値 | トレンド・季節成分 |
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つの不確かな情報を合わせると、どちらか一方より確かな値が得られる点が要です。
次の時刻には、この20.8℃(分散0.8)から再び予測します。時間が経つ間に温度が変わる可能性があるので、予測の分散には「モデルの不確かさ」Qが足されて0.8より大きくなります。そこへ次の観測が来て、また分散の比で混ぜます。この「予測で分散が増え、観測で縮む」の繰り返しがカルマンフィルタの全体です。観測の分散が予測より大きい場合はKが0.5未満になり、推定値は予測寄りに留まります。
予測ステップと更新ステップの式:状態空間モデルの行列版
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の点検に使います。予測の不確かさPが大きいほど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点ですでに観測そのままより悪化しました。カルマンフィルタは速度も状態として持ち、「今はこのあたりまで進んでいるはず」と予測してから混ぜるので、この実験では後ろ向き移動平均より小さい誤差で現在位置を推定できました。ただし、モデルが外れると追従遅れも生じます。移動平均との選択は、許容する遅れとモデル化の可否に応じて判断します。
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の過小設定による失敗:急な速度変化への追従遅れ
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年9月時点の各公式ページで確認した値です。
| ライブラリ | 最新版(公開日) | 任せられること |
|---|---|---|
| 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です。UnobservedComponents などの状態空間モデルは、カルマンフィルタで尤度を計算してノイズの分散を最尤推定します。MATLABでは、Control System Toolboxの kalman 関数がプラントモデルとQ・Rからフィルタを設計します。制御系の開発環境としての全体像はMATLAB/Simulinkの解説にまとめています。
カルマンフィルタが向かない場面と外れ値への備え
カルマンフィルタは前提が崩れると、何も警告を出さずに間違った値を返します。採用を見送るか、手当てが要るのは次の場面です。
- 外れ値が混ざる観測:ガウス雑音を仮定しているため、1回の大きな外れ値に推定値が引っ張られます。後述のゲートで対処します。
- 状態の変化をモデルにできない:状態変化を妥当なモデルで表せない場合、推定精度を保証できません。ゲインは設定した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で一定のスカラー量を直接観測し、観測誤差が無相関かつ等分散なら、事前情報を持たない極限で最小二乗解(標本平均)と一致します。状態が動く場合に広げたものがカルマンフィルタです。
カルマンフィルタの状態方程式とは何ですか?
状態が1ステップでどう変わるかを表す式 x_k = F x_(k−1) + w_k です。Fが変化の規則、w_kがモデルで表しきれないゆらぎで、その共分散がQです。観測との関係は観測方程式 z_k = H x_k + v_k で別に書きます。
カルマンフィルタと拡張カルマンフィルタの違いは何ですか?
カルマンフィルタは状態の変化と観測が線形の場合の手法で、拡張カルマンフィルタ(EKF)は非線形の式を現在の推定値のまわりで線形化して同じ計算を使います。EKFはApollo計画の航法研究の中で生まれました。
カルマンゲインとは何ですか?
観測と予測の差を状態の修正量へ変換する係数です。状態を同じ単位で直接観測する1次元モデル(H=1)では 予測の分散 / (予測の分散 + 観測の分散) で、0〜1の値を取ります。1に近いほど観測を、0に近いほど予測を信じます。
Pythonでカルマンフィルタを使うにはどうすればよいですか?
線形ならNumPyで数十行で書けます。ライブラリならfilterpy(最新版は2018年の1.4.5)、pykalman、OpenCVの cv2.KalmanFilter、ノイズ分散も推定したい時系列ならstatsmodelsが選択肢です。