lazy diary

統計とその周辺

ローカルレベルモデルとカルマンフィルタ

ちょっとだけ話題(?)になっていたカルマンフィルタを参考文献の本で学んでいます。話題になっているから学んでいるわけではないですが、時系列分析の基本だと思うので。

今回は第2章のローカルレベルモデルに対するカルマンフィルタを実装します。ローカルレベルモデルとは次のような最もシンプルな状態空間モデルのことです。
観測値をy_t、観測できない潜在変数を\alpha_tとして
f:id:mt19:20190519101423j:plain
という仮定を置いたモデルです。ただし、これだと初期値について何も言及がないので
f:id:mt19:20190519101430j:plain
のように初期分布を置いてあげます。従って未知パラメータは\sigma^2_{\epsilon}, \sigma^2_{\eta}, a_1, P_1の4つです。未知パラメータのうち初期分布に関わるa_1, P_1については、何らかの事前知識があれば良いのですが今回は無いのでいわゆる無情報事前分布にしてしまいます。つまり、分散P_1を大きくしてしまいます。このときa_1はもはやどこでも変わらないので、常識的な範囲で適当に置いてもらって構いません。他のパラメータについては最尤推定で決める方法があるようですが、今回は既知として自分で設定しました。

使うデータは参考文献の本から提供されている体重の推移データです。データを使う場合は出版社のページから飛んでください。
カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― / 野村 俊一 著 | 共立出版

ローカルレベルモデルはシンプルであり、stanを使えばあっという間に実装できそうな気がしますが、今回はあくまでカルマンフィルタのお勉強です。ということで、アルゴリズムを導出して、その式を使って実装しています。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

#load data and plot
weight_data = pd.read_csv("./Weight.dat", header = None)

plt.figure(figsize = (8, 4))
plt.plot(weight_data)
plt.xlabel("day")
plt.ylabel("weight(kg)")
plt.grid(True)
plt.show()

f:id:mt19:20190519103005p:plain

ローカルレベルモデルのカルマンフィルタを実装です。

class LocalLevelModel():
    """
    this class gives LocalLevelModel using KalmanFilter.
    model accepts only univariate pandas data.
    if missing data is described np.NaN, model deals with it automatically
    """
    def __init__(self, data, a_ini = 0.0, P_ini = 1000, ep_ini = 1.0, eta_ini = 0.1):
        self.data = data.values
        self.n = len(data)
        self.a_ini = a_ini
        self.P_ini = P_ini
        self.ep_ini = ep_ini
        self.eta_ini = eta_ini
        self.missing_ind = np.where(np.isnan(self.data) == True)[0]
        self.dropna_data = self.data[~np.isnan(self.data)]
    
    
    #calculate filtered estimator
    def fit(self):
        v, F, a_con, a = np.empty(self.n), np.empty(self.n), np.empty(self.n), np.empty(self.n)
        P, P_con, K, L = np.empty(self.n), np.empty(self.n), np.empty(self.n), np.empty(self.n)

        #初期値代入
        a[0] = self.a_ini
        P[0] = self.P_ini
        
        
        #逐次計算
        for t in range(self.n - 1):
            F[t] = P[t] + self.ep_ini ** 2
            K[t] = P[t] / F[t]
            L[t] = 1 - K[t]
            
            if t in self.missing_ind:
                a_con[t] = a[t]
                P_con[t] = P[t]
                a[t + 1] = a_con[t]
                P[t + 1] = P_con[t] + self.eta_ini ** 2
            else:
                v[t] = self.data[t] - a[t]
                a_con[t] = a[t] + K[t] * v[t]
                P_con[t] = P[t] * L[t]
                a[t + 1] = a_con[t]
                P[t + 1] = P_con[t] + self.eta_ini ** 2
        
        t = t + 1
        F[t] = P[t] + self.ep_ini ** 2
        K[t] = P[t] / F[t]
        L[t] = 1 - K[t]
        if t in self.missing_ind:
            a_con[t] = a[t]
            P_con[t] = P[t]
        else:
            v[t] = self.data[t] - a[t]
            a_con[t] = a[t] + K[t] * v[t]
            P_con[t] = P[t] * L[t]
        
        return a_con, P_con, F, v, L
    
    #state smoothing
    def smoothing_fit(self):
        #smoothed state
        hat_alpha = np.empty(self.n)
        #smoothed state variance
        V = np.empty(self.n)
        
        r = np.empty(self.n)
        N = np.empty(self.n)
        
        #initialize
        r[self.n - 1] = 0
        N[self.n - 1] = 0
        
        a_con, P_con, F, v, L = self.fit()
        
        #逐次計算
        for t in range(self.n - 1):
            if t in self.missing_ind:
                r[self.n - 2 - t] = r[self.n - 1 - t]
                hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
                
                N[self.n - 2 - t] = N[self.n - 1 - t]
                V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
            else:
                r[self.n - 2 - t] = v[self.n - 1 - t] / F[self.n - 1 - t] + L[self.n - 1 - t] * r[self.n - 1 - t]
                hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
            
                N[self.n - 2 - t] = 1 / F[self.n - 1 - t] + (L[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
                V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
            
        t = t + 1
        hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
        V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
        
        return hat_alpha, V

変数名は参考文献に寄っているので、本が無いと分かりにくいかもしれません。重要な変数については説明を置いているつもりです。欠測値(np.NaN)にも自動で対応するようにしました。

fit()ではフィルタ化推定量とその分散、smoothing_fit()では平滑化状態とその分散を返してくれます。

まずは体重の推移データで、欠測値が一切無いときをプロットしてみます。
f:id:mt19:20190526120150p:plain
f:id:mt19:20190526120158p:plain

平滑化状態は全ての観測データを利用した予測なので、フィルタ化推定量に比べて幅が狭くなっているのが分かります。

次に欠測値を混ぜてみます。20日から40日の間は体重を測定しなかったと仮定して、データをモデルに投げます。
f:id:mt19:20190526120227p:plain
f:id:mt19:20190526120234p:plain
欠測値のあたりだけ不自然に膨らんでいます。データが無いので幅が広くなっていることが分かります。

上の欠測値のあるときのフィルタ化推定量というのは実質n期先の予測に他なりません。ローカルレベルモデルでは直前の状態を参考にするしかないため、そうすると「直前の値をひたすらまっすぐ取り続けるであろう」というなんとも情けない予測になってしまいます。ここが良くないですね。

おまけ1(2019/05/25)

ためしに外れ値を混ぜてみたらどういう結果になるか調べてみる。
f:id:mt19:20190526120315p:plain
f:id:mt19:20190526120323p:plain
上がフィルタ化推定量、下が平滑化状態の推定量です。思いのほか引っ張られていない?よく分からぬ。

おまけ2(2019/05/25)

TFPでフィルタ化推定量を求めてみた。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import tensorflow_probability as tfp
import tensorflow as tf
import tensorflow.contrib.eager as tfe
tfe.enable_eager_execution()
    
tfl = tf.linalg
tfd = tfp.distributions


weight_data = pd.read_csv("./Weight.dat", header = None)
weight_data = weight_data.values

num_dims = 1
num_timesteps = weight_data.shape[0]
step_mu = 0.0
step_std = 1.0
noise_mu = 0.0
noise_std = 1.0
ini_loc = tf.convert_to_tensor(np.array([0], dtype = np.float32))
ini_scale = tf.convert_to_tensor(np.array([1000]), dtype = np.float32)
scale = tf.convert_to_tensor(np.array([0.1]), dtype = np.float32)
model = tfd.LinearGaussianStateSpaceModel(
    num_timesteps = num_timesteps,
    transition_matrix = tfl.LinearOperatorIdentity(num_dims) ,
    transition_noise = tfd.MultivariateNormalDiag(loc = tf.zeros([num_dims]), scale_diag = scale),
    observation_matrix = tfl.LinearOperatorIdentity(num_dims),
    observation_noise = tfd.MultivariateNormalDiag(loc = tf.zeros([num_dims]), scale_diag = tf.ones([num_dims])),
    initial_state_prior = tfd.MultivariateNormalDiag(loc = ini_loc, scale_diag = ini_scale)
)

input_data = tf.convert_to_tensor(weight_data.astype(np.float32))
output_data = model.forward_filter(input_data)
prediction = np.array(output_data[1]).flatten()
std = np.sqrt(np.array(output_data[2]).flatten())
x = range(len(weight_data))
z = stats.norm.ppf(q = 0.975)
plt.figure(figsize = (8, 4))
plt.scatter(x, weight_data)
plt.plot(weight_data)
plt.plot(prediction)
plt.fill_between(x, prediction - z * std, prediction + z * std, alpha = 0.2)
plt.grid(True)
plt.show()

f:id:mt19:20190526121640p:plain
たぶん同じように動いてくれているはず・・・。
もしコピーしてjupyterなどで自分で動かしてみようという物好きな方がいれば注意してください。kernelを再起動してまっさらな状態でeager_executionを実行しないとエラーが出ます。
python - tf.enable_eager_execution must be called at program startup ONLY in SPYDER IDE - Stack Overflow

参考文献

カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)

この記事については第2章を参考にしました。