[PYTHON] Parameter estimation with Kalman filter

In the weather forecast, the Kalman filter is used to estimate the parameters to convert the model output into forecast variables (for example, maximum temperature). The method is described in the Japan Meteorological Agency (1997).

Parameter estimation

When there is no time dependence

Here, we try to estimate $ a $ and $ b $ in the linear equation $ y = ax + b $. Use the $ y_i = a x_i + b + v_i $ formula to generate simulated data and examine it.

Let $ a = 2 $, $ b = 5 $, $ x $ be a uniform random number from -5 to 5, and $ v_i $ give a Gaussian random number with mean 0 and standard deviation 2.

KLM.py

def update(P, C, R, x_hat, obs, I):
  """
  P:Error covariance matrix
  C:Observation system coefficient matrix
  R:Observed noise variance matrix
  """
  #Kalman gain
  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)

To illustrate the estimated $ a $ and $ b $, Then, in a short time, it can be seen that $ a $ and $ b $ are estimated.

If there is time dependence

When the model changes on the way or when the parameters change seasonally. This time, when 2 is added to each coefficient on the 180th day.

KLM2.py

import numpy as np

def update(P, C, R, x_hat, obs, I):
  """
  P:Error covariance matrix
  C:Observation system coefficient matrix
  R:Observed noise variance matrix
  """
  #Kalman gain
  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)

It followed the change brilliantly.

References

Japan Meteorological Agency (1997) Latest Numerical Weather Prediction Guide Chapter 5 Application System 1 Kalman Filter pp.66 --78