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 $, KLM.png 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)

KLM2.png

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

Recommended Posts

Parameter estimation with Kalman filter
Change point detection with Kalman filter
Non-linear Kalman filter (1)
Parameter tuning with luigi (2)
Parameter tuning with luigi
Create filter with scipy
Use DATE_FORMAT with SQLAlchemy filter
[Python] Filter spreadsheets with gspread
Kalman filter that you can understand
HMM parameter estimation implementation in python
Make a filter with a django template
Motor parameter estimation from corresponding points