PDF 形式のマニュアルはダウンロードページからダウンロードできます。

本書は、「Monte Carlo Dynamic Classifier ツール開発および実験支援作業」において作成したプログラムの実行手順を説明するものです。 Monte Carlo Dynamic Classifier ツールは、任意の観測データ系列のモデル推定、およびその推定モデルによる状態系列の推定を行うプログラムです。 推定されたモデルは、異なる観測データ系列に適用してモデルの尤度を計算することにより、観測データ系列のクラス分類などに応用できます。

MCDC ツールは、以下のプログラム群から構成されます。

MCDCTrain

モデル推定プログラム

MCDCTest

モデルの尤度計算プログラム

Graphs

推定されたモデルのグラフ描画関数群

本書の以下の各章では、これらのプログラムの実行方法と実行例、プログラム構成 について説明します。
また、本ツールに付属するサンプルプログラムの一部では、 CMU Graphics Lab Motion Capture Database にて公開されているモーションキャプチャデータを利用します。 このデータの利用許諾条件についてライセンスに記載します。

MCDC ツールでは、観測データ系列に対してモデル推定を行う MCDCTrain と、推定モデルを用いて未知の観測データ系列の状態系列の推定を行う MCDCTest が提供されます。 このほかに、推定モデルをグラフ描画する関数群を集めたツールとして Graphs クラスが提供されます。本章ではこれらのプログラムの実行方法について説明します。

MCDCTrain 関数は、観測データ系列に対してモデル推定を行うプログラムです。以下のように実行します。

[ IDX, SKP, OKP, FV, GV, XE, YE, loglik ] = MCDCTrain( ...
  algorithm, ...
  grids, ...
  stateKernelGens, ...
  obsKernelGens, ...
  stateMeanFuncs, ...
  obsMeanFuncs, ...
  gridDimForGM, ...
  splineHandle, ...
  x0, ...
  xaux, ...
  u, ...
  y', ...
  N, ...
  J, ...
  K, ...
  aspect ...
);

また、以前の実行で出力されたファイルを読み込み、反復実行を継続することもできます。この場合は MCDCTrain 関数を以下のように実行します。 継続実行の際には、指定された matfile からパラメータを読み込み、前回と同じ設定で実行を再開します。ただし、 Name, Value の組を指定することで、前回の設定を上書きして実行できます。

[IDX, SKP, OKP, FV, GV, XE, YE, loglik] = MCDCTrain(matfile, aspect, Name, Value. ...)

MCDCTrain 関数のパラメータは以下のとおりです。

パラメータ	データ型	内容
algorithm	Algorithm	アルゴリズム種別
xGrids	G double[] cell	状態空間上に構成される格子点の座標。各次元の値を double 配列として格納したセル配列として表す
stateKernelGens	Dx handle cell	状態遷移関数のカーネル事前分布。 Dx は状態空間のうち推定対象とされる次元数とする
obsKernelGens	P handle cell	観測関数のカーネル事前分布
gridDimForGramMatrix	int	グラム行列作成に用いる次元
splineHandle	handle	スプライン補間に用いるアルゴリズム。現在は GenericSpline のハンドルに固定
x0	Dx*1 double	状態変数の初期状態
xaux	Da*T double	追加状態データ
u	D*T double	制御データ
z	P*T double	観測データ
N	int	Particle filter 実行時の粒子数
J	int	MCMC iteration の全体の反復回数
K	int	観測オフセット
aspect	Aspect	システム設定情報 (ログ出力先等)

algorithm パラメータには、以下のいずれかのクラスを選択できます。

クラス名	アルゴリズム内容
PSMC	状態遷移関数と観測関数をガウス過程からのランダムサンプリングで生成する。パラメータ推定は行わない
PMCMC	状態遷移関数と観測関数をガウス過程からのランダムサンプリングで生成する。パラメータ推定に Particle marginal Metropolis-Hastings 法を用いる
PMCMC2	状態遷移関数と観測関数をガウス過程からのランダムサンプリングで生成する。パラメータ推定に Particle marginal Metropolis-Hastings 法を用いる。状態空間の特定の一次元のみを共分散関数に用いる

アルゴリズムとして PMCMC2 を利用する場合には、さらに、平均値関数と共分散関数の学習方式、共分散関数に用いるカーネルパラメータの学習方式を以下から選択できます。

クラス名	アルゴリズム内容
MCDCStrategyChoice1	平均値関数を学習せず、常にアンカーモデルを利用する。共分散関数にカーネル関数を用いる
MCDCStrategyChoice2	現在の GP surface を平均値関数として関数をサンプリングする。共分散関数にカーネル関数を用いる
MCDCStrategyChoice3	現在の GP surface を平均値関数として関数をサンプリングする。共分散関数にカーネル関数を用いず、固定された共分散行列を用いる
RBFKernelGeneratorStrategyChoice1	RBF カーネルのカーネルパラメータを事前分布からランダムサンプリングする
RBFKernelGeneratorStrategyChoice2	RBF カーネルのカーネルパラメータを現在値からのランダムウォークによって生成する

MCDCTrain 関数の実行結果として戻される値は、以下のとおりです。

値	データ型	内容
IDX	J*1 double	j 回目の反復までの累積受理回数
SKP	A*D*2 double	第 a 受理反復 における、状態遷移関数のカーネルパラメータ (sigma, l)
OKP	A*P*2 double	第 a 受理反復における、観測関数のカーネルパラメータ (sigma, l)
FV	A*D*G double	第 a 受理反復における、状態遷移関数による格子点上の値。 G は格子点のサイズに一致する多次元配列
GV	A*P*G double	第 a 受理反復における、観測関数による格子点上の値。 G は格子点のサイズに一致する多次元配列
XE	A*T*D double	第 a 受理反復における、時刻 t の状態変数の推定平均
YE	A*T*P double	第 a 受理反復における、時刻 t の観測変数の推定平均
loglik	A*1 double	第 a 受理反復における、 particle filter による推定の対数尤度

実行時には、 aspect の設定にしたがって中間状態が保存されます。これは .mat 形式のダンプファイルとして出力されます。ファイルに格納されているデータは以下のとおりです。

値	データ型	内容
algorithm	Algorithm	MCDCTrain のパラメータに渡された algorithm
in	MCDCInput	MCDCTrain のパラメータのうち、 algorithm と aspect を除くすべて
out	MCDCOutput	MCDCTrain の戻り値のすべて (途中状態) および、現在の反復回数 j

実行中には、 aspect の設定にしたがってログファイルが出力されます。これはテキスト形式のファイルになります。ログファイルは以下のような形式になります。

2014/04/26 11:29:36 - Iteration 18 / 200
2014/04/26 11:29:36 -   StateKernel[1]: [ Sigma=8.932992, L=7.409991 ]
2014/04/26 11:29:36 -   StateKernel[2]: [ Sigma=4.410584, L=2.304895 ]
2014/04/26 11:29:36 -   ObsKernel[1]: [ Sigma=1.711191, L=9.248074 ]
2014/04/26 11:29:36 -   ObsKernel[2]: [ Sigma=9.867704, L=9.599577 ]
2014/04/26 11:29:36 -   ObsKernel[3]: [ Sigma=7.700470, L=8.679293 ]
2014/04/26 11:29:36 -   Creating GramMatrix using 1 dim...
2014/04/26 11:29:36 -     StateKernel[1] Done
2014/04/26 11:29:36 -     StateKernel[2] Done
2014/04/26 11:29:36 -     ObsKernel[1] Done
2014/04/26 11:29:36 -     ObsKernel[2] Done
2014/04/26 11:29:36 -     ObsKernel[3] Done
2014/04/26 11:29:36 -   Drawing GP surface...
2014/04/26 11:29:36 -   Estimating using particle filter... (N=500)
2014/04/26 11:30:20 -   Acceptance log probability = 232745.510426
2014/04/26 11:30:20 -   logLH = -5347068.232758, accepted
2014/04/26 11:30:20 - Elapsed time is 44.324939 seconds.
2014/04/26 11:30:20 -   j: 8 bytes
2014/04/26 11:30:20 -   IDX: 1600 bytes
2014/04/26 11:30:20 -   SKP: 128 bytes
2014/04/26 11:30:20 -   OKP: 192 bytes
2014/04/26 11:30:20 -   FV: 61504 bytes
2014/04/26 11:30:20 -   GV: 92256 bytes
2014/04/26 11:30:20 -   XE: 122752 bytes
2014/04/26 11:30:20 -   YE: 184128 bytes
2014/04/26 11:30:20 -   loglik: 32 bytes

MCDCTest 関数は、 MCDCTrain によって得られたモデルを用いて未知のデータの状態を推定します。複数の異なるモデルを用いて状態推定を行い、それぞれの尤度を比較することで、クラス分 類問題への応用も可能です。以下のように実行します。

[ result, FnState, FnObs ] = MCDCTest(u, y, N, modelFile)

MCDCTest 関数のパラメータは以下のとおりです。

パラメータ	データ型	内容
u	D*T double	制御データ
y	P*T double	観測データ
N	int	Particle filter 実行時の粒子数
modelFile	chars	モデルファイル

MCDCTest 関数の実行結果として戻される値は、以下のとおりです。

値	データ型	内容
result	ParticleFilter	j 回目の反復までの累積受理回数
FnState	handle	状態遷移関数
FnObs	handle	観測関数

戻り値の result に、 particle filter による状態推定の結果が含まれます。これは以下の構造を持ちます。

値	データ型	内容
Particles	N*T*D double	時刻 t における各粒子の座標
Weights	N*T double	時刻 t における各粒子の重み
Loglik	double	推定された状態の対数尤度

MCDCTrain によって推定されたモデルは、 Graph クラスに定義された関数群を用いて PDF ファイルに出力することができます。以下のグラフを出力可能です。

• Graphs.YE
• Graphs.XE
• Graphs.YEMean
• Graphs.XEMean
• Graphs.Loglik
• Graphs.Rmse

観測データ系列と、推定モデルを用いた粒子フィルタによって観測データを追跡した結果の時間推移グラフを出力します。 MCMC iteration の中で、指定された特定の反復における推定モデルを用いた結果を出力します。観測データ系列の次元ごとにグラフが出力されます。 以下のように実行します。

Graphs.YE( ...
  outputFileNamePrefix, ...
  matFileName, ...
  iterations, ...
  times ...
);

パラメータは以下のとおりです。

パラメータ	データ型	内容
outputFileNamePrefix	char[]	出力グラフファイルのファイル名接頭辞。ファイル名は、指定された接頭辞に、次元数と拡張子 .pdf を付け加えたものになります
matFileName	char[]	推定モデルを含む Matlab データファイル名
iterations	int[]	プロットする推定値の MCMC iteration を列挙した配列
times	int[]	プロットするデータ系列の範囲。省略時はデータ系列の全体を出力します

推定モデルを用いた粒子フィルタによる状態系列推定値の時間推移グラフを出力します。 MCMC iteration の中で、指定された特定の反復における推定モデルを用いた結果を出力します。状態系列の次元ごとにグラフが出力されます。 以下のように実行します。

Graphs.XE( ...
  outputFileNamePrefix, ...
  matFileName, ...
  iterations, ...
  times ...
);

パラメータは以下のとおりです。

パラメータ	データ型	内容
outputFileNamePrefix	char[]	出力グラフファイルのファイル名接頭辞。ファイル名は、指定された接頭辞に、次元数と拡張子 .pdf を付け加えたものになります
matFileName	char[]	推定モデルを含む Matlab データファイル名
iterations	int[]	プロットする推定値の MCMC iteration を列挙した配列
times	int[]	プロットするデータ系列の範囲。省略時はデータ系列の全体を出力します

2.3.3. Graphs.YEMean

Graphs.YEMean( ...
  outputFileNamePrefix, ...
  matFileName, ...
  times ...
);

Graphs.XEMean( ...
  outputFileNamePrefix, ...
  matFileName, ...
  times ...
);

Graphs.Loglik( ...
  outputFileNamePrefix, ...
  matFileName1, ...
  matFileName2, ...
  ...
);

Graphs.Rmse( ...
  outputFileNamePrefix, ...
  matFileName ...
);

[x, y] = KitagawaModel(1000, 0.5, 28, 8, 0.6, 30, 10, 0.05, 0.06, 0.07, 0.08, 0.1, 0.1);
u  = repmat(cos(1.2 * [1:T]), 2, 1)';

grids = { ...
  [-30:2:30], ...
  [-30:2:30]  ...
};

stateKernelGens = { ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10))  ...
};

obsKernelGens = {
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10)), ...
  RBFKernelGenerator(UniformDistribution(0.01, 10), UniformDistribution(0.01, 10))  ...
};

stateMeanFuncs = { ...
  @(x) (a1 .* x(1,:) + b1 .* x(1,:) / (1 + x(1,:) .^ 2)), ...
  @(x) (a2 .* x(2,:) + b2 .* x(2,:) / (1 + x(2,:) .^ 2))  ...
};

obsMeanFuncs = { ...
  @(x) (d1 .* x(1,:) .^ 2 + d2 .* x(2,:) .^ 2), ...
  @(x) (d3 .* x(1,:) .^ 2                    ), ...
  @(x) (                    d4 .* x(2,:) .^ 2)  ...
};

kernelGeneratorStrategy = @RBFKernelGeneratorStrategyChoice2;
mcdcStrategy = @MCDCStrategyChoice2;
algorithm = PMCMC2(kernelGeneratorStrategy, mcdcStrategy);
gridDimForGramMatrix = 1;
splineHandle = @GenericSpline;

x0 = zeros(1, 2);
N  = 500;
J  = 100000;

aspect = Aspect();
aspect.LogFileName = 'logs/KitagawaModelPMCMC2.log';
aspect.MatFileNamePrefix = 'logs/KitagawaModelPMCMC2';
aspect.SavesIntermediateMat = true;
aspect.IntermediateMatInterval = 100;

[ IDX, SKP, OKP, FV, GV, XE, YE, loglik ] = MCDCTrain( ...
  algorithm, …
  grids, ...
  stateKernelGens, ...
  obsKernelGens, ...
  stateMeanFuncs, ...
  obsMeanFuncs, ...
  gridDimForGramMatrix, ...
  splineHandle, ...
  x0, ...
  [], ... % No auxiliary states
  u', ...
  y', ...
  N, ...
  J, ...
  0, ...
  aspect ...
);