「カルマンフィルターのパワーを暴露する」
「カルマンフィルターのパワーを大公開！」
データサイエンティストとして、将来の値を予測するためにトレンドをモデル化する必要があることがあります。統計的または機械学習ベースのアルゴリズムに焦点を当てる誘惑がありますが、私は異なるオプションを提案するためここにいます：カルマンフィルタ（KF）です。
1960年代初頭、ルドルフ・E・カルマンはKFで複雑なシステムをモデル化する方法を革新しました。航空機や宇宙船を目的地に案内したり、スマートフォンがこの世界で場所を見つけるのを可能にすることから、このアルゴリズムはデータと数学を組み合わせて信じられないほどの正確さで将来の状態の推定値を提供します。
このブログでは、カルマンフィルタがどのように動作するかについて詳しく解説し、Pythonの例を示し、この技術の真の力を強調します。単純な2Dの例から始めて、コードを変更してより高度な4Dの空間に適応する方法、そして洗練された前身である拡張カルマンフィルタをカバーする最後の部分まで進みます。予測アルゴリズムとフィルタの世界への旅に参加してください。
カルマンフィルタの基本
KFは、時間にわたって収集された観測と他の測定からのノイズと過去の状態の統計的分布を表す共分散行列のセットを構築し、継続的に更新することで、システムの状態の推定を提供します。他の既存のアルゴリズムとは異なり、システムと外部ソースとの数学的関係を定義することでソリューションを直接拡張および改良することが可能です。このプロセスはかなり複雑で入り組んでいるように聞こえるかもしれませんが、手順は予測と更新に要約することができます。これらのフェーズは連携して、システムの状態推定を反復的に修正および改善します。
予測ステップ：
このフェーズでは、モデルの既知の事後推定および時間Δkのステップに基づいて、システムの次の将来の状態を予測することに関連しています。数学的には、状態空間の推定値を次のように表します：
ここで、Fは状態推移行列であり、制御入力とプロセスノイズに関係なく、状態が1つのステップから別のステップに進化する方法をモデル化します。行列Bは、制御入力uₖが状態に与える影響をモデル化します。
次の状態の推定値とともに、アルゴリズムは共分散行列Pによって表される推定値の不確実性も計算します：
予測された状態の共分散行列は、予測の信頼性と精度を表し、システム自体からのプロセスノイズの共分散行列Qの影響を受けます。この行列を更新ステップの後続の方程式に適用して、カルマンフィルタがシステムに保持している情報を補正し、将来の状態推定を改善します。
更新ステップ：
更新ステップでは、カルマンゲイン、状態推定値、および共分散行列を更新します。カルマンゲインは、新しい測定が状態推定値にどれだけの影響を与えるかを決定します。計算には、状態と受け取る予定の測定値との関係を示す観測モデル行列Hと、測定値の誤差の測定ノイズ共分散行列Rが含まれます：
要するに、Kは予測の不確実性を測定値の中に存在する不確実性とバランスさせようとします。前述のように、カルマンゲインは状態推定値と共分散を補正するために適用されます。それぞれの方程式は次のように示されます：
括弧内の計算は、状態推定の残差であり、実際の測定値とモデルが予測した測定値との差です。
カルマンフィルタの真の美しさは、新しい情報が受け取られると状態と共分散の両方を連続的に更新する再帰的な性質にあります。これにより、ノイズのある観測値のストリームを受けるシステムを統計的に最適な方法で時間の経過とともに洗練することができます。
The Kalman filter in action
カルマンフィルタの根底にある方程式には非常に圧倒されることがあり、数学のみでその動作を完全に理解するには状態空間の理解が必要です（このブログの範囲外です）。しかし、いくつかのPythonの例を使って生き生きと説明を試みます。最も単純な形式では、カルマンフィルタオブジェクトを次のように定義できます：
import numpy as npclass KalmanFilter: """ 線形動的システムのクラシックなカルマンフィルタの実装 カルマンフィルタはノイズのある観測値からシステムの状態を推定するための最適な再帰型データ処理アルゴリズムです。 属性： F（np.ndarray）：状態遷移行列。 B（np.ndarray）：制御入力行列。 H（np.ndarray）：観測行列。 u（np.ndarray）：制御入力。 Q（np.ndarray）：プロセスノイズの誤差共分散行列。 R（np.ndarray）：測定ノイズの誤差共分散行列。 x（np.ndarray）：前のステップ（k-1）の平均状態推定値。 P（np.ndarray）：前のステップ（k-1）の状態誤差共分散。 """ def __init__(self, F, B, u, H, Q, R, x0, P0): """ 必要な行列と初期状態を持つカルマンフィルタを初期化します。 パラメータ： F（np.ndarray）：状態遷移行列。 B（np.ndarray）：制御入力行列。 H（np.ndarray）：観測行列。 u（np.ndarray）：制御入力。 Q（np.ndarray）：プロセスノイズの誤差共分散行列。 R（np.ndarray）：測定ノイズの誤差共分散行列。 x0（np.ndarray）：初期状態推定。 P0（np.ndarray）：初期状態誤差共分散行列。 """ self.F = F # 状態遷移行列 self.B = B # 制御入力行列 self.u = u # 制御ベクトル self.H = H # 観測行列 self.Q = Q # プロセスノイズの誤差共分散 self.R = R # 測定ノイズの誤差共分散 self.x = x0 # 初期状態推定 self.P = P0 # 初期推定誤差共分散 def predict(self): """ 次のステップの状態と状態誤差の予測を行います。 """ self.x = self.F @ self.x + self.B @ self.u self.P = self.F @ self.P @ self.F.T + self.Q return self.x def update(self, z): """ 最新の測定値で状態推定を更新します。 パラメータ： z（np.ndarray）：現在のステップの測定値。 """ y = z - self.H @ self.x S = self.H @ self.P @ self.H.T + self.R K = self.P @ self.H.T @ np.linalg.inv(S) self.x = self.x + K @ y I = np.eye(self.P.shape[0]) self.P = (I - K @ self.H) @ self.P return self.x非線形システムの課題上記のpredict()関数とupdate()関数を使用して、前述の手順を反復します。上記のフィルタの設計は、任意の時系列で機能し、推定が実際の値とどのように比較されるかを示すため、次の単純な例を構築しましょう：
import numpy as npimport matplotlib.pyplot as plt# 再現性のためにランダムシードを設定np.random.seed(42)# オブジェクトの実際の位置のグラウンドトゥルースをシミュレートtrue_velocity = 0.5 # 単位ごとの時間ステップnum_steps = 50time_steps = np.linspace(0, num_steps, num_steps)true_positions = true_velocity * time_steps# ノイズを含んだ測定値をシミュレートmeasurement_noise = 10 # 測定値をノイズ化するためにこの値を増やすnoisy_measurements = true_positions + np.random.normal(0, measurement_noise, num_steps)# 真の位置とノイズのある測定値をプロットplt.figure(figsize=(10, 6))plt.plot(time_steps, true_positions, label='真の位置', color='green')plt.scatter(time_steps, noisy_measurements, label='ノイズのある測定値', color='red', marker='o')plt.xlabel('時間ステップ')plt.ylabel('位置')plt.title('真の位置とノイズのある測定値の時間変化')plt.legend()plt.show()
実際のところ、「真の位置」は不明ですが、ここでは参照のためにプロットします。「ノイズのある測定値」は、カルマンフィルタに供給される観測点です。非常に基本的な行列のインスタンス化を行い、カルマンモデルはカルマンゲインの適用により迅速に収束するため、ある程度は問題ありませんが、特定の状況ではモデルへのウォームスタートを実行することが合理的かもしれません。
# Kalman Filterの初期化F = np.array([[1, 1], [0, 1]]) # 状態遷移行列B = np.array([[0], [0]]) # 制御入力なしu = np.array([[0]]) # 制御入力なしH = np.array([[1, 0]]) # 観測関数Q = np.array([[1, 0], [0, 3]]) # プロセスノイズの共分散R = np.array([[measurement_noise**2]]) # 観測ノイズの共分散x0 = np.array([[0], [0]]) # 初期の状態推定値P0 = np.array([[1, 0], [0, 1]]) # 初期の推定誤差共分散kf = KalmanFilter(F, B, u, H, Q, R, x0, P0)# 推定位置と速度の保存用の配列を作成estimated_positions = np.zeros(num_steps)estimated_velocities = np.zeros(num_steps)# Kalman Filterのループfor t in range(num_steps): # 予測 kf.predict() # 更新 measurement = np.array([[noisy_measurements[t]]]) kf.update(measurement) # フィルタリングされた位置と速度を保存 estimated_positions[t] = kf.x[0] estimated_velocities[t] = kf.x[1]# 真の位置、ノイズのある観測値、およびKalman Filterの推定値をプロットplt.figure(figsize=(10, 6))plt.plot(time_steps, true_positions, label='真の位置', color='green')plt.scatter(time_steps, noisy_measurements, label='ノイズのある観測値', color='red', marker='o')plt.plot(time_steps, estimated_positions, label='Kalman Filterの推定値', color='blue')plt.xlabel('タイムステップ')plt.ylabel('位置')plt.title('時間にわたるKalman Filterの推定値')plt.legend()plt.show()
この非常にシンプルなソリューションの設計でも、モデルはノイズを貫通して真の位置を見つけるので、まあまあの結果が得られます。これは単純なアプリケーションでもうまく動作しますが、トレンドはより微妙で、外部イベントの影響を受けることがしばしばあります。これを処理するためには、状態空間表現だけでなく、新しい情報が到着したときに推定値を計算し、共分散行列を修正する方法を変更する必要があります。別の例でこれを詳しく見てみましょう。
4次元での移動物体の追跡
空間と時間で物体の移動を追跡したいとします。さらに、角度の回転による力が加わるというより現実的な例をシミュレートします。このアルゴリズムが高次元の状態空間表現に対してどれだけ適応性があるかを示すために、線形力を仮定しますが、実際にはそうではありません（この後により現実的な例を探求します）。以下のコードは、この特定のシナリオに対してKalman Filterを修正する方法の例です。
import numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dclass KalmanFilter: """ 線形ダイナミックシステムのクラシックなカルマンフィルターの実装。 カルマンフィルターは、ノイズのある観測からシステムの状態を推定するための最適な再帰型データ処理アルゴリズムです。 Attributes: F (np.ndarray): 状態遷移行列。 B (np.ndarray): 制御入力行列。 H (np.ndarray): 観測行列。 u (np.ndarray): 制御入力。 Q (np.ndarray): プロセスノイズの共分散行列。 R (np.ndarray): 観測ノイズの共分散行列。 x (np.ndarray): 前のステップ（k-1）の状態の平均推定値。 P (np.ndarray): 前のステップ（k-1）の状態の推定誤差共分散行列。 """ def __init__(self, F=None, B=None, u=None, H=None, Q=None, R=None, x0=None, P0=None): """ 必要な行列と初期状態でカルマンフィルターを初期化します。 Parameters: F (np.ndarray): 状態遷移行列。 B (np.ndarray): 制御入力行列。 H (np.ndarray): 観測行列。 u (np.ndarray): 制御入力。 Q (np.ndarray): プロセスノイズの共分散行列。 R (np.ndarray): 観測ノイズの共分散行列。 x0 (np.ndarray): 初期状態推定値。 P0 (np.ndarray): 初期推定誤差共分散行列。 """ self.F = F # 状態遷移行列 self.B = B # 制御入力行列 self.u = u # 制御入力 self.H = H # 観測行列 self.Q = Q # プロセスノイズの共分散行列 self.R = R # 観測ノイズの共分散行列 self.x = x0 # 初期状態推定 self.P = P0 # 初期推定誤差共分散 def predict(self): """ 次のタイムステップの状態と状態の共分散を予測します。 """ self.x = np.dot(self.F, self.x) + np.dot(self.B, self.u) self.P = np.dot(np.dot(self.F, self.P), self.F.T) + self.Q def update(self, z): """ 最新の観測値で状態推定を更新します。 Parameters: z (np.ndarray): 現在のステップの観測値。 """ y = z - np.dot(self.H, self.x) S = np.dot(self.H, np.dot(self.P, self.H.T)) + self.R K = np.dot(np.dot(self.P, self.H.T), np.linalg.inv(S)) self.x = self.x + np.dot(K, y) self.P = self.P - np.dot(np.dot(K, self.H), self.P)# シミュレーション用パラメータtrue_angular_velocity = 0.1 # タイムステップごとのラジアン半径 = 20num_steps = 100dt = 1 # タイムステップ# タイムステップの作成time_steps = np.arange(0, num_steps*dt, dt)# 真の状態true_x_positions =
ここで注目すべきいくつかの興味深い点があります。上のグラフでは、観測の反復が始まると、モデルが推定された真の状態に迅速に修正される様子がわかります。モデルはまた、システムの真の状態を適切に特定する能力もあります。推定値が真の軌道と交差する（「真の軌道」）ことが示されています。この設計は、一部の現実世界のアプリケーションに適しているかもしれませんが、システムに非線形な力が作用する場合には適していません。代わりに、カルマンフィルタの別のアプリケーションである拡張カルマンフィルタを考慮する必要があります。拡張カルマンフィルタは、到着観測の非線形性を線形化するために、ヤコビアン行列を解に組み込み、テイラー級数展開を実行することによって、最後とは異なるアルゴリズムです。この拡張を数学的に表現するために、主なアルゴリズムの計算は次のようになります:
状態予測用の式、fは前の状態推定に適用される非線形状態遷移関数で、uは前のタイムステップの制御入力です。
誤差共分散予測の式、Fは状態遷移関数のヤコビアンであり、Pは前の誤差共分散であり、Qはプロセスノイズの共分散行列です。
k時点での計測値zの観測式、hは状態予測に適用される非線形観測関数であり、vは観測ノイズの加算です。
カルマンゲインの更新計算式、Hは観測関数のヤコビアンであり、Rは測定ノイズの共分散行列です。
非線形観測関数とカルマンゲインを組み込んだ状態推定の修正計算、最後に誤差共分散を更新する式:
最後の例では、これにより特定のオブジェクトの角度回転の非線形影響を線形化するためにヤコビアンを使用し、コードを適切に修正します。拡張カルマンフィルタの設計は、線形近似の初等的な仮定が状態推定に誤りをもたらす可能性があるため、カルマンフィルタよりも複雑です。さらに、ヤコビアンの計算は複雑になり、計算負荷が高く、特定の状況を定義するのが困難である場合もあります。ただし、正しく設計された場合、拡張カルマンフィルタの実装は通常、カルマンフィルタを上回る性能を発揮します。
前回のPythonの例に基づいて、EKFの実装を以下に示します:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
class ExtendedKalmanFilter:
"""
The Extended Kalman Filter (EKF) implementation.
This filter is suitable for systems with non-linear dynamics by linearising
the system model at each time step using the Jacobian.
Attributes:
state_transition (callable): The state transition function for the system.
jacobian_F (callable): Function to compute the Jacobian of the state transition.
H (np.ndarray): The observation matrix.
jacobian_H (callable): Function to compute the Jacobian of the observation model.
Q (np.ndarray): The process noise covariance matrix.
R (np.ndarray):
シンプルな要約
このブログでは、カルマンフィルタ（KF）の構築と適用方法、および拡張カルマンフィルタ（EKF）の実装方法について詳しく説明しました。さて、これらのモデルの使用例、利点、欠点について要約しましょう。
KF：このモデルは、状態の遷移と観測行列が状態の線形関数（一部ガウスノイズを伴う）であると仮定できる線形システムに適用できます。このアルゴリズムは、次のような場合に適用することが考えられます：
- 一定速度で移動するオブジェクトの位置と速度を追跡する
- ノイズが確率的であるか線形モデルで表現できる信号処理アプリケーション
- 線形モデルでモデル化できる経済予測
KFの主な利点は、（行列計算に従えば）アルゴリズムがかなり単純であり、他のアプローチよりも計算量が少なく、時間内のシステムの真の状態の非常に正確な予測と推定を提供できることです。欠点は、典型的には複雑な現実世界のシナリオでは線形性の仮定が成り立たないことです。
EKF：実質的には、EKFをKFの非線形版とみなすことができ、ヤコビアンの使用によってサポートされています。EKFは以下のような場合に考慮されます：
- 測定値とシステムダイナミクスが通常非線形であるロボティクスシステム
- 非一定速度または角速度の変化がある追跡航空機や宇宙船などを追跡するトラッキングおよびナビゲーションシステム。
- 最も近代的な「スマート」な車に見られるクルーズコントロールやレーンキーピングの実装時の自動車システム。
EKFは、非線形システムに対してKFよりも優れた推定値を提供することが多いですが、ヤコビアン行列の計算により計算量が増加する可能性があります。このアプローチはまた、テイラー級数展開からの1次線形近似を前提としていますが、高度に非線形なシステムに対しては適切な仮定とは言えない場合があります。ヤコビアンの追加により、モデルの設計がより困難になる場合もあり、利点にもかかわらず、シンプルさと相互運用性のためにKFを実装することがより適切かもしれません。
特に記載がない限り、すべての画像は著者によるものです。
We will continue to update VoAGI; if you have any questions or suggestions, please contact us!
Was this article helpful?
93 out of 132 found this helpful
Related articles
