はじめに

前回に引き続きカルマンフィルターの勉強を少しだけ進めました。
今回もこちらを参考にさせていただき、同じ事を Python でやってみたいと思います。
ローカルレベルモデルのパラメータの最尤推定となります。

例によって、初学者が勉強しながら残す備忘録ですので、誤り等含む可能性があることご了承ください。同じレベルから同じ事を勉強する人の助けになればという位置づけです。

ローカルレベルモデル

事前準備

まず準備として下記のパッケージを import します。

import matplotlib.pyplot as plt
import scipy.stats
import scipy.optimize
import numpy as np
import math

統計関連はRの方が充実しているようですが、今回の内容は scipy があれば十分です。
パッケージがインストールされていない場合は

pip install numpy
pip install scipy
pip install matplotlib

などのコマンドでインストールできます。

シミュレーションデータの生成

ローカルレベルモデルで扱うデータは、いわゆるランダムウォークするデータだと理解しています。真値は1ステップ毎に正規分布誤差に基づいて変動していきます。

また、観測においては真値をそのまま取得することはできず、こちらも正規分布に基づく観測誤差が加わった結果しか得られないものとします。

変動誤差と観測誤差には相関は存在しない(互いに独立したノイズである)前提で次のプログラムで示すような生成モデルとしています。

このようなモデルが裏に隠れている前提で、観測値のみからモデルのパラメータを推定していくことで、最も尤もらしい予測値を得ることを目標とします(モデル自体は予測できている前提になります)。

def create_random_walk(n, v, w):
    """ランダムウォークの真値と観測データを生成
    
    ランダムウォークの真値と観測誤差入り観測データを生成

    Args:
        n(int): 生成するデータ数
        v(float): 変動誤差の分散
        w(float): 観測誤差の分散
    
    Returns:
        List[flot], List[flot]: 真値,観測データ
    """
    x = np.cumsum(np.random.normal(0.0, math.sqrt(v), n))
    y = x + np.random.normal(0.0, math.sqrt(w), n)
    return x, y

NumPy にある cumsum() という関数を今回はじめて知りましたが、その要素より前の要素を積算した値を格納してくれます。離散系での積分を行うのにとても便利そうです。

ローカルレベルモデルでの予測

下記が今回のメインとなるローカルレベルモデルの予測部分です。
当初、カルマンゲインの求め方がなぜこの式になるのか私にはさっぱりわからなかったのですが、その点について初心者なりに式を分解してみたのが前回の記事となりますので、よろしければご参照ください。

def predict_llm(x, x_var, y, v, w):
    """ローカルレベルモデルでの予測
    
    ローカルレベルモデルで計測値から次のデータを予測

    Args:
        x(float): 現在の予測値
        x_var(float): 現在の予測値の尤度の分散
        y(float): 新しい観測地
        v(float): 変動誤差の分散
        w(float): 観測誤差の分散

    Returns:
        flot, float, float: 新しい予測値、新しい予測値の尤度の分散値、尤度
    """
    
    # 尤度
    likelihood = scipy.stats.norm.pdf(y, loc=x, scale=np.sqrt(x_var + v + w))
    
    # 現在値の予測
    predict_x     = x           # 期待値は前回の値のそのまま
    predict_x_var = x_var + v   # 分散を加算される
    
    # カルマンゲインを計算
    kalman_gain = predict_x_var / (predict_x_var + w)
    
    # 計測値と合成
    next_x     = predict_x + kalman_gain * (y - predict_x)
    next_x_var = (1 - kalman_gain) * predict_x_var
    
    # 2変数の最適合成という点では上はこう書いても等価
#   next_x     = (w * predict_x +  predict_x_var * y) / (predict_x_var + w)
#   next_x_var = (predict_x_var * w) / (predict_x_var + w)
    
    return next_x, next_x_var, likelihood

ここで likelihood を求めている箇所が非常に重要です。
１ステップ分の計算ごとに、現在の予測値とその確からしさを求めていくわけですが、前回予測した x の値と x のバラツキ予想であるx_var に、さらに今回追加されたであろう変動誤差 v と観測誤差 w を加えて、得られた確率密度関数(正規分布と仮定)に対して、新しい観測値yが得られる確率を求めています。
このように発生した事象から逆方向にそれが起こる尤もらしさを求計算したのが尤度と呼ばれるものです。
そして、その尤度を最大化するパラメータを求めるのが最尤推定となります。

なお、ここで分散値を単純に足し算できるところが偏差ではなく分散を使うとうれしい事象の１つです。
誤差を持った要因が重ねあわされると、誤差が互いに独立である限りにおいて、分散値を加算した誤差が加わったのと同じことになります。

(余談：誤差を組み合わせると分散が加算されるというのはとても重要な法則です。ちなみにこれはホワイトノイズの乗った同じものを観測した信号を100個足し合わせるとS/Nがその平方根をとって偏差に戻すことで10倍の改善になるるという信号処理などでも利用するおなじみの法則です。)

最尤推定を行ってみる

シミュレーションデータ生成

下記のようにシミュレーションデータを作成し、これを真値とすることにします。

v_true = 2.0  # 変動誤差の分散(後で推定したい定数1)
w_true = 6.0  # 観測誤差の分散(後で推定したい定数2)
x_true, y = create_random_walk(1000, v_true, w_true)

計算関数の準備

まず下記のように、任意の誤差分散 v, w を行い推定を行う関数を準備します。

def calc_predict(y, v, w):
    """ ローカルレベルモデルで計測値から予測を実施
    """
    n = len(y)
    x = np.ndarray(n)
    x_var = np.ndarray(n)
    likelihood = np.ndarray(n)
    x[0] = 0
    x_var[0] = 0
    likelihood[0] = 1.0
    for i in range(1, n):
        x[i], x_var[i], likelihood[i] = predict_llm(x[i-1], x_var[i-1], y[i], v, w)
    return x, likelihood

また、結果の評価用にRMSEを計算する関数も作って見ます。
これは私の興味本位なのでこの指標が妥当かは不明です(すみません)。

def calc_rmse(x_true, y, v, w):
    """ 予測値と真値のRMSEを計算してみる 
    """
    x, _ = calc_predict(y, v, w)
    return np.sqrt(np.mean((x - x_true)**2))

ためしに誤差の分散が機知の前提で計算

ためしに推定に先立って、生成したデータを誤差の分散が機知の前提で計算してみます。

x, likelihood = calc_predict(y, v_true, w_true)
print('RMSE = %f' % calc_rmse(x_true, y, v_true, w_true))

# 最後の100個を表示
plt.figure(figsize=(12, 3), dpi=100)
plt.title('[v=%f w=%f] RMSE=%f' % (v_true, w_true, rmse))
plt.plot(x_true[9900:10000], c='red', label='ground truth')
plt.plot(y[9900:10000], c='green', label='observation')
plt.plot(x[9900:10000], c='blue', label='prediction')
plt.legend()
plt.savefig('plot.png')
plt.show()

こんな感じの結果になりました。

誤差パラメータの最尤推定

早速ですが、scipy.optimize.minimize を使って、尤度が最大になるパラメータを予想してみます。
最大になる値を求めるのにminimizeを使うので、評価関数の戻り値を符号反転しています。

def minimize_likelihood(x):
    """ scipy.optimize の関数
        分散は負の値を取らないので exp して使う
    """
    _, likelihood = calc_predict(y, np.exp(x[0]), np.exp(x[1]))
    return -np.sum(np.log(likelihood))  # 尤度最大化したいので符号を反転させる

# 予測実施
result = scipy.optimize.minimize(minimize_likelihood, x0=[1.0, 1.0], method="L-BFGS-B")

v_predict = np.exp(result.x[0])
w_predict = np.exp(result.x[1])
print('v_predict = %f (v_true = %f)' % (v_predict, v_true))
print('w_predict = %f (w_true = %f)' % (w_predict, w_true))
print('RMSE = %f' % calc_rmse(x_true, y, v_predict, w_predict))

結果

v_predict = 2.055061 (v_true = 2.000000)
w_predict = 6.005451 (w_true = 6.000000)
RMSE = 1.612539

ちゃんとモデルのパラメータが予測できているように思います。

尤度は確率ですので、各ステップで起こりうる尤度をすべて掛け合わせたものが全体の尤度になります。
ここですべてのステップの掛け算を行うと値が大きくなりすぎるので、対数を取って扱うのが一般的だそうです。前人の方々の知恵はすばらしいですね。
対数を取ると乗算は加算に変わりますので総和を取るようにしています。
また、分散値は正の値しか取りません。scipy.optimizeには制約を掛ける機能もあるようですが、先人に習って、expをとることにします。
ここは-∞～+∞の値が単調増加性を保ったまま0～+∞に写像できる関数であれば何でもよいのかとは思います。

真値が既知の前提でRMSE最小化の計算をしてみる(おまけ)

ここからは私の興味本位ですが、真値がわかっている前提でRMSE(Root Mean Squared Error)が最小になるように最適化を行ってみます。

def minimize_mse(x):
    """ scipy.optimize の関数
    """
    return calc_rmse(x_true, y, np.exp(x[0]), np.exp(x[1]))

# 予測の実施
result = scipy.optimize.minimize(minimize_mse, x0=[1.0, 1.0], method="L-BFGS-B")

v_predict = np.exp(result.x[0])
w_predict = np.exp(result.x[1])
print('v_predict = %f (v_true = %f)' % (v_predict, v_true))
print('w_predict = %f (w_true = %f)' % (w_predict, w_true))
print('RMSE = %f' % calc_rmse(x_true, y, v_predict, w_predict))

結果

v_predict = 1.558741 (v_true = 2.000000)
w_predict = 4.740400 (w_true = 6.000000)
RMSE = 1.612352

なんと、推定値の真値に対する二乗誤差自体はわずかに減少したものの元のパラメータの推定誤差は非常に大きくなっています。

目的は本来のモデルのパラメータを推定することですので、ここでも最尤推定の強力さを垣間見た気がします。

#おわりに

この手の分野の専門ではないので、時々思い出したように挑戦しますが、なかなか勉強も進んでいません。
しかし真値がわかっている前提であっても素人アプローチではモデルの推定をするのが難しいというのは私としては新しい発見でした。
また進展あれば記事を書いてみたいと思います。

なお、今回のコードは、こちらに置いています。