Bei Wettervorhersagen werden Parameter mithilfe des Kalman-Filters geschätzt, um die Modellausgabe in Vorhersagevariablen (z. B. Maximaltemperatur) umzuwandeln. Die Methode ist in der Meteorological Agency (1997) beschrieben.
Hier schätzen wir $ a $ und $ b $ in der linearen Gleichung $ y = ax + b $. Verwenden Sie die Formel $ y_i = a x_i + b + v_i $, um simulierte Daten zu generieren und zu untersuchen.
Sei $ a = 2 $, $ b = 5 $, $ x $ eine einheitliche Zufallszahl von -5 bis 5, und $ v_i $ ergibt eine Gaußsche Zufallszahl mit einem Mittelwert von 0 und einer Standardabweichung von 2.
KLM.py
def update(P, C, R, x_hat, obs, I):
"""
P:Fehlerkovarianzmatrix
C:Zahlenmatrix des Beobachtungssystems
R:Beobachtungsrauschdispersionsmatrix
"""
#Kalman gewinnen
G = P * C / (C.T * P * C + R)
x_hat = x_hat + G * (obs - C.T * x_hat)
P = (I - G * C.T) * P
return x_hat, P
a = 2.
b = 5.
x = np.random.uniform(-5, 5, 365)
v = np.random.normal(0, 2, 365)
y = []
for i in range(365):
y.append(a * x[i] + b + v[i])
A = np.mat([1])
P = np.mat([[1, 0], [0, 1]])
R = np.mat([1])
I = np.identity(2)
x_hat = np.mat([[0], [0]])
for i in range(365):
C = np.mat([[x[i]], [1]])
obs = np.mat([y[i]])
x_hat, P = update(P, C, R, x_hat, obs, I)
Um die geschätzten $ a $ und $ b $ zu veranschaulichen, Es ist ersichtlich, dass $ a $ und $ b $ in kurzer Zeit geschätzt werden.
Wenn sich das Modell unterwegs ändert oder wenn sich die Parameter saisonal ändern. Diesmal, wenn am 180. Tag zu jedem Koeffizienten 2 addiert wird.
KLM2.py
import numpy as np
def update(P, C, R, x_hat, obs, I):
"""
P:Fehlerkovarianzmatrix
C:Zahlenmatrix des Beobachtungssystems
R:Beobachtungsrauschdispersionsmatrix
"""
#Kalman gewinnen
G = P * C / (C.T * P * C + R)
x_hat = x_hat + G * (obs - C.T * x_hat)
P = (I - G * C.T) * P
return x_hat, P
a = 2.
b = 5.
x = np.random.uniform(-5, 5, 365)
v = np.random.normal(0, 2, 365)
y = []
for i in range(365):
if i < 180:
y.append(a * x[i] + b + v[i])
else:
y.append((a + 2) * x[i] + (b + 2) + v[i])
A = np.mat([1])
P = np.mat([[1, 0], [0, 1]])
R = np.mat([1])
S = 0.1
I = np.identity(2)
x_hat = np.mat([[0], [0]])
for i in range(365):
C = np.mat([[x[i]], [1]])
obs = np.mat([y[i]])
P = P + S ** 2
x_hat, P = update(P, C, R, x_hat, obs, I)
Es folgte der Veränderung brillant.
Meteorological Agency (1997) Neuester Leitfaden zur numerischen Vorhersage Kapitel 5 Anwendungssystem 1 Kalman-Filter S. 66 - 78
Recommended Posts