[PYTHON] Kalman Filter, den Sie verstehen können
Über diesen Artikel
Lassen Sie uns den Kalman-Filter am Beispiel des Random Walks implementieren.
Was ist ein Kalman-Filter?
Es ist eine Technik, um die Person mit dem Namen Noise zu entfernen und das wahre Erscheinungsbild abzuschätzen.
Theoretisch
Betrachten Sie als typische Einstellung ein System wie das folgende.
\begin{align}
x_{t+1} &= F_t x_t + G_t w_t \tag{1.1}\\
y_t &= H_t x_t + v_t \tag{1.2}\\
\end{align}
$ x_t $ repräsentiert den internen Zustand des Systems und $ y_t $ ist seine Beobachtung. Der interne Zustand $ x_t $ entwickelt sich im Laufe der Zeit gemäß $ (1.1) $. Andererseits wird der beobachtete Wert $ y_t $ gemäß $ (1.2) $ mit Rauschen überstrichen. Der Zweck besteht darin, aus dem beobachteten Wert $ y_t $ auf einen internen Zustand $ x_t $ zu schließen, der nicht direkt beobachtet werden kann. Nehmen wir der Einfachheit halber an, dass $ w_t $ und $ v_t $ unabhängigen Gaußschen Verteilungen folgen und nur gleichzeitig eine Autokorrelation aufweisen. Verwendung der Unabhängigkeit von $ w_t $ und $ v_t $ und des Bayes'schen Theorems
\begin{align}
p(x_t|Y^t) &= \frac{p(y_t|x_t)p(x_t|Y^{t-1})}{p(y_t|Y^{t-1})} \tag{2.1}\\
p(x_{t+1}|Y^t) &= \int p(x_{t+1}|x_t)p(x_t|Y^t)dx_t \tag{2.2}
\end{align}
Der relationale Ausdruck wird erhalten. Wobei $ Y ^ t $ der Wert von $ y_t $ bis zum Zeitpunkt $ t $ ist
Y^t=\left\{y_1, y_2, ..., y_t\right\}
Nehme an, dass
Jetzt,
Andererseits, wenn die beobachteten Werte bis zur Endzeit erhalten werden
p(x_t|Y^N) = \int \frac{p(x_{t+1}|x_t)p(x_t|Y^t)p(x_{t+1}|Y^N)}{p(x_{t+1}|Y^t)}dx_{t+1} \tag{2.3}
Kann auch abgeleitet werden. Hier,
Experimentieren Sie mit Random Walk
Zielloser Spaziergang
Betrachten Sie als einfaches Beispiel den Fall eines zufälligen Spaziergangs. Die Formeln für $ (1.1) $ und $ (1.2) $ sind
\begin{align}
x_{t+1} &= x_t + w_t \tag{3.1}\\
y_t &= x_t + v_t \tag{3.2}
\end{align}
Es wird sein. Es wird jedoch angenommen, dass $ w_t $ und $ v_t $ unabhängigen Gaußschen Verteilungen folgen und nur gleichzeitig eine Autokorrelation aufweisen. Sie können sich $ (3.1) $ als die wahre Bewegung des Zufallslaufs und $ (3.2) $ als das zusätzliche Rauschen vorstellen, um die beobachteten Werte zu erhalten. Ich kenne nur das beobachtete $ y_t $, aber ich möchte das ursprüngliche $ x_t $ erraten, indem ich das Rauschen entferne.
Übrigens, machen wir einen zufälligen Spaziergang mit Python. Die folgenden drei Pakete sind ausreichend.
import numpy as npimport numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context('poster')
Wir werden die obige Formel implementieren. Machen Sie zuerst einen zufälligen Spaziergang von $ (3.1) $ und fügen Sie dann Rauschen wie $ (3.2) $ hinzu.
def random_walker(start_position=0, mean=0, deviation=1, n_steps=99, seed=None):
if seed is not None:
np.random.seed(seed=seed)
move = np.random.normal(loc=mean, scale=deviation, size=n_steps)
position = np.insert(move, 0, start_position)
position = np.cumsum(position)
return position
def add_noise(position, mean=0, deviation=10, seed=None):
if seed is not None:
np.random.seed(seed=seed)
n_observation = len(position)
noise = np.random.normal(loc=mean, scale=deviation, size=n_observation)
observation = position + noise
return observation
Lassen Sie uns tatsächlich einen lauten zufälligen Spaziergang erzeugen. Die Startposition ist 0. Außerdem sei die tatsächliche Random-Walk-Dispersion ($ w_t $ Dispersion) 1 und die Rauschdispersion ($ v_t $ Dispersion) 10. Der Durchschnitt ist 0 für beide. Stellen Sie den Zeitschritt auf 1 ~ 100 ein.
true_position = random_walker(start_position=0, mean=0, deviation=1, n_steps=99, seed=0)
observed_position = add_noise(true_position, mean=0, deviation=10, seed=0)
Lassen Sie uns ein Diagramm erstellen.
plt.plot(true_position, 'r--', label='True Positions')
plt.plot(observed_position, 'y', label='Observed Ppositions')
plt.title('Random Walk')
plt.xlabel('time step')
plt.ylabel('position')
plt.legend(loc='best')
Es ist ein wunderbarer zufälliger Spaziergang. True Positions ist der wahre Wert des Random Walks $ x_t $, und Observed Positions ist die verrauschte Beobachtung $ y_t $.
Der beobachtete Wert ist verschmutzt.
Kalman Filter
Schreiben wir die Kalman-Filterformel für dieses Random-Walk-Modell auf. Angenommen, $ w_t $ folgt einer Gaußschen Verteilung mit einem Mittelwert von 0 und einer Varianz von $ q $, und $ v_t $ folgt einer Gaußschen Verteilung mit einem Mittelwert von 0 und einer Varianz von $ r
\begin{align}
&\hat{x}_{t/t} = \hat{x}_{t/t-1} + K_t(y_t - \hat{x}_{t/t-1}) \\
&\hat{x}_{t+1/t} = \hat{x}_{t/t} \\
&K_t = \frac{P_{t/t-1}}{P_{t/t-1} + r} \\
&P_{t/t} = \frac{r P_{t/t-1}}{P_{t/t-1} + r} \\
&P_{t+1/t} = P_{t/t} + q \\
\end{align}
Die Formel zum Finden von $ \ hat {x} \ _ {t | N} $ in umgekehrter Reihenfolge der Zeit lautet
\begin{align}
&\hat{x}_{t/N} = \hat{x}_{t/t} + C_t(\hat{x}_{t+1/N} - \hat{x}_{t+1/t}) \\
&C_t = \frac{P_{t/t}}{P_{t/t} + q} \\
&P_{t/N} = P_{t/t} + C^2_t(P_{t+1/N} - P_{t+1/N})
\end{align}
Es wird sein. Lassen Sie uns dies implementieren.
class Simple_Kalman:
def __init__(self, observation, start_position, start_deviation, deviation_true, deviation_noise):
self.obs = observation
self.n_obs = len(observation)
self.start_pos = start_position
self.start_dev = start_deviation
self.dev_q = deviation_true
self.dev_r = deviation_noise
self._fit()
def _forward(self):
self.x_prev_ = [self.start_pos]
self.P_prev_ = [self.start_dev]
self.K_ = [self.P_prev_[0] / (self.P_prev_[0] + self.dev_r)]
self.P_ = [self.dev_r * self.P_prev_[0] / (self.P_prev_[0] + self.dev_r)]
self.x_ = [self.x_prev_[0] + self.K_[0] * (self.obs[0] - self.x_prev_[0])]
for t in range(1, self.n_obs):
self.x_prev_.append(self.x_[t-1])
self.P_prev_.append(self.P_[t-1] + self.dev_q)
self.K_.append(self.P_prev_[t] / (self.P_prev_[t] + self.dev_r))
self.x_.append(self.x_prev_[t] + self.K_[t] * (self.obs[t] - self.x_prev_[t]))
self.P_.append(self.dev_r * self.P_prev_[t] / (self.P_prev_[t] + self.dev_r))
def _backward(self):
self.x_all_ = [self.x_[-1]]
self.P_all_ = [self.P_[-1]]
self.C_ = [self.P_[-1] / (self.P_[-1] + self.dev_q)]
for t in range(2, self.n_obs + 1):
self.C_.append(self.P_[-t] / (self.P_[-t] + self.dev_q))
self.x_all_.append(self.x_[-t] + self.C_[-1] * (self.x_all_[-1] - self.x_prev_[-t+1]))
self.P_all_.append(self.P_[-t] + (self.C_[-1]**2) * (self.P_all_[-1] - self.P_prev_[-t+1]))
self.C_.reverse()
self.x_all_.reverse()
self.P_all_.reverse()
def _fit(self):
self._forward()
self._backward()
Wenden wir es auf den zuvor erwähnten zufälligen Spaziergang an. In der Praxis kennen Sie normalerweise nicht den wahren Wert von zufälligen Spaziergängen oder Geräuschstreuung und müssen ihn abschätzen. Lassen Sie uns hier zunächst den wahren Wert verwenden.
kf = Simple_Kalman(observed_position, start_position=0, start_deviation=1, deviation_true=1, deviation_noise=10)
Lassen Sie uns ein Diagramm erstellen.
plt.plot(true_position, 'r--', label='True Positions')
plt.plot(observed_position, 'y', label='Observed Ppositions')
plt.plot(kf.x_, 'blue' ,label='Foward Estimation')
plt.plot(kf.x_all_, 'black', label='Smoothed Estimation')
plt.title('Random Walk')
plt.xlabel('time step')
plt.ylabel('position')
plt.legend(loc='best')
Vorwärtsoptimierung
Es ist schön geworden.
Verweise
[Nichtlinearer Kalman-Filter](https://www.amazon.co.jp/%E9%9D%9E%E7%B7%9A%E5%BD%A2%E3%82%AB%E3%83%AB%E3% 83% 9E% E3% 83% B3% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF-% E7% 89% 87% E5% B1% B1-% E5% BE% B9 / dp / 4254201486 /) Asakura Buchhandlung, Toru Katayama
Die Formeln basieren auf diesem Buch. Es ist eine wundervolle Abdeckung mit entferntem Lärm.
Recommended Posts