[PYTHON] Parameterschätzung mit Kalman-Filter

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.

Parameter Schätzung

Wenn es keine Zeitabhängigkeit gibt

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 es Zeitabhängigkeit gibt

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.

Verweise

Meteorological Agency (1997) Neuester Leitfaden zur numerischen Vorhersage Kapitel 5 Anwendungssystem 1 Kalman-Filter S. 66 - 78