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).
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.
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.
Japan Meteorological Agency (1997) Latest Numerical Weather Prediction Guide Chapter 5 Application System 1 Kalman Filter pp.66 --78
Recommended Posts