Monte Carlo Dynamic Classifier

本ページでは、MCDC ツールによるモデル推定の実行例として、Kitagawa モデルに従う観測データ系列のモデル推定を行う手順を紹介します。

Kitagawa モデルは、非線形状態空間モデルの一つで、時刻 t に対して以下の式で表されます。

状態遷移関数

x1(t) = a1 * x1(t-1) + b1 * x1(t-1) / (1 + x1(t-1) ^ 2) + c1 * cos(1.2 * t) + sqrt(q1) * randn;
x2(t) = a2 * x2(t-1) + b2 * x2(t-1) / (1 + x2(t-1) ^ 2) + c2 * cos(1.2 * t) + sqrt(q2) * randn;

観測関数

y1(t) = d1 * x1(t) ^ 2 + d2 * x2(t) ^ 2 + sqrt(r1) * randn;
y2(t) = d3 * x1(t) ^ 2                  + sqrt(r2) * randn;
y3(t) =                  d4 * x2(t)     + sqrt(r3) * randn;

ただし a, b, c, d はモデルパラメータであり、q, r は系列に加えられるノイズの大きさを制御するパラメータです。

この実験では、Kitagawa モデルから生成された観測データ系列を与えて、モデル推定を行い、推定されたモデルの観測データ系列への当てはめを確認します。モデルパラメータは以下のように設定しました。

実験に用いたパラメータは以下のとおりです。

@(x) (a1 .* x(1,:) + b1 .* x(1,:) / (1 + x(1,:) .^ 2)), ...
@(x) (a2 .* x(2,:) + b2 .* x(2,:) / (1 + x(2,:) .^ 2))  ...

@(x) (d1 .* x(1,:) .^ 2 + d2 .* x(2,:) .^ 2), ...
@(x) (d3 .* x(1,:) .^ 2                    ), ...
@(x) (                    d4 .* x(2,:)     )  ...

推定モデルを用いて観測データ系列を推定した結果を以下のグラフに示します。Kitagawa モデルでは時刻ごとに観測値が大きく変化するため、グラフでは時刻 t = 950 から 1000 の範囲を拡大して示しています。赤色の実線が観測データ、青色の実線が推定モデルによる推定平均です。青色の破線は推定平均±標準偏差の幅を表します。

観測データ系列の推定結果 (順に y1, y2, y3)

MCMC iteration ごとの推定モデルの対数尤度を次のグラフに示します。

推定モデルの対数尤度

この実験を行った計算機環境と実行時間は以下のとおりです。

Kitagawa model の推定では、真のモデルから cos(1.2 * t) の項を除いたものをアンカーモデルとしています。 別ページの Lorenz model の結果と比較すると、観測データ系列に対する当てはまりは悪く、上記の項の影響をうまく扱えていないことが考えられます。