[PYTHON] Estimation des paramètres avec le filtre de Kalman

Dans les prévisions météorologiques, l'estimation des paramètres est effectuée à l'aide du filtre de Kalman pour convertir la sortie du modèle en variables de prévision (par exemple, la température maximale). La méthode est décrite dans l'Agence météorologique (1997).

Estimation des paramètres

Quand il n'y a pas de dépendance au temps

Ici, estimons $ a $ et $ b $ dans l'équation linéaire $ y = ax + b $. Utilisez la formule $ y_i = a x_i + b + v_i $ pour générer des données simulées et les examiner.

Soit $ a = 2 $, $ b = 5 $, $ x $ un nombre aléatoire uniforme de -5 à 5, et $ v_i $ donne un nombre aléatoire gaussien avec une moyenne de 0 et un écart type de 2.

KLM.py

def update(P, C, R, x_hat, obs, I):
  """
  P:Matrice de covariance d'erreur
  C:Matrice des numéros du système d'observation
  R:Matrice de dispersion du bruit observé
  """
  #Gain de Kalman
  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)

Pour illustrer les estimés $ a $ et $ b $, On peut voir que $ a $ et $ b $ sont estimés sur une courte période de temps.

S'il y a une dépendance au temps

Lorsque le modèle change en cours de route ou lorsque les paramètres varient selon les saisons. Cette fois, lorsque 2 est ajouté à chaque coefficient le 180e jour.

KLM2.py

import numpy as np

def update(P, C, R, x_hat, obs, I):
  """
  P:Matrice de covariance d'erreur
  C:Matrice des numéros du système d'observation
  R:Matrice de dispersion du bruit observé
  """
  #Gain de Kalman
  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)

Il a suivi le changement avec brio.

Les références

Agence météorologique (1997) Dernier guide des prévisions numériques Chapitre 5 Système d'application 1 Filtre de Kalman pp.66 --78