Monte Carlo Dynamic Classifier

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

Lorenz 方程式は非線形方程式の一つで、ここでは、Lorenz 方程式をもとに以下の非線形状態空間モデルを考えます。

状態遷移関数

x1(t) = (Sigma * x2(t-1) * dt + x1(t-1) * (1 - Sigma * dt))   + sqrt(q1) * randn;
x2(t) = (x1(t-1) * (Rho - y1(t-1))) * dt + x2(t-1) * (1 - dt) + sqrt(q2) * randn;

観測関数

y1(t) = (x1(t-1) * x2(t-1) * dt) + y1(t-1) * (1 - Beta * dt)  + sqrt(r1) * randn;
y2(t) = Alpha0 * x1(t-1) .^ 2                                 + sqrt(r2) * randn;
y3(t) = Alpha0 * x2(t-1) .^ 2                                 + sqrt(r3) * randn;

ただし Sigma, Rho, Beta, Alpha0 はモデルパラメータであり、q, r は系列に加えられるノイズの大きさを制御するパラメータです。dt は時間間隔です。

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

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

@(x) (Sigma .* x(2,:) .* dt + x(1,:) .* (1 - Sigma .* dt)),  ...
@(x) (x(1,:) .* (Rho - x(3,:)) .* dt + x(2,:) .* (1 - dt))   ...

@(x) ((x(1,:) .* x(2,:) .* dt) + x(3,:) .* (1 - Beta * dt)), ...
@(x) (Alpha0 * x(1,:) .^ 2),                                 ...
@(x) (Alpha0 * x(2,:) .^ 2)                                  ...

推定モデルを用いて観測データ系列を推定した結果を以下のグラフに示します。赤色の実線が観測データ、青色の実線が推定モデルによる推定平均です。青色の破線は推定平均±標準偏差の幅を表します。

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

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

推定モデルの対数尤度

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

Lorenz model の推定では、真のモデルに等しいアンカーモデルを設定しているため、2 回目以降の MCMC iteration の結果はすべて棄却されています。 推定モデルは観測データ系列を比較的よく表現できていますが、これは真のモデルをアンカーモデルにしているためと考えられます。