という訳で、混合正規分布についてEMアルゴリズムのサンプルプログラムを作ってみましょうか。
で、どうせなら、ちょっとメンドクサイですが、多変量の場合についてのものに取り組んでみましょう。
なお、あくまでも「サンプル」で、実用に耐え得るものかは分かりません（汗
ちょっと試して遊ぶためだけのものを目指しましょう。
でも、EMアルゴリズムで遊ぶ人ってどんな人なんでしょう・・・
さて、ここまでにあげたEMアルゴリズムの更新式を見ての通り、プログラム中のあちこちで和をとったりすることになるのですが、その際のインデックスの使い方を決めておきましょう。
基本、更新式の時の指標の使い方に準ずるとして
１．データ番号のインデックスを示すのには int i を使う。また、データの総数は size_t ndata とする。
２．ベクトルの成分を示すのには int j を使う。また、ベクトルの次元は size_t ndim とする。
３．クラスタの番号を示すのには int k を使う。また、クラスタの総数は size_t ncls とする。
クラスタ、とは今の場合混合正規分布を構成する各正規分布の事です。
ともすると（書いてる人が）忘れがちですが、もともとクラスタリングでブログを書いてたはずなので（汗）、以降クラスタという言葉を使いましょう。
また、size_t ndata、 ndim、 ncls を別な構造体に保持させる場合なども、その構造体のメンバ名を ndata、 ndim、 ncls として対応付けやすいようにしましょう。
さて、多変量の混合正規分布ですが、前にも載せたように密度関数は
となります。
ここで x は M 次元のベクトルです。
で、ここまでの書き方に従えば M = ndim です。
一変量ならパラメータは平均、分散のスカラー値2個だけですが、多変量の場合
μ : 平均ベクトル : ndim 次元ベクトル
Σ : 分散共分散行列 : ndim 行 x ndim 列の正方行列
となります。
さらに密度関数を計算するには、
|Σ| : 行列Σのデターミナント
Σ^-1 : 行列Σの逆行列
が必要になります。
密度関数だけでもなんか大がかりになりますね・・・
で、これらのベクトルや行列は、全て1次元配列で表現しましょう。
行列なんかは2次元配列の方が分かり易いと思うのですが、そうすると1変量のときに困ったことになりそうなのです。。。
但し行列は row major で。
BLASなんかは col major ですが、あれは大嫌いです。。。
つまり
という行列は、 {0, 1, 2, 3} というように、列方向に順にデータを入れてゆきます。
混合分布の場合、正規分布をクラスタ数 ncls 個分考えなければなりません。
なので上のパラメータが ncls 個分必要になる訳ですが、そのためにまず、正規分布1個のパラメータを格納する構造体を用意しましょう。
その上で、混合分布を表す構造体を別途用意し、それに正規分布パラメータの構造体を ncls 個持たせるようにしましょう。
正規分布パラメータの構造体メンバとしては
その正規分布の
平均ベクトル
分散共分散行列
分散共分散行列の行列式
分散共分散行列の逆行列
でしょうか。
またこれらを1次元配列で表すなら要素数も持たせなければなりませんね。
平均ベクトルの次元も、分散共分散行列の行数、列数も全てデータベクトルの次元 ndim に一致するので、構造体にもこの次元数も持たせておきます。
以上から、多変量正規分布のパラメータを格納する構造体として以下を用意します。
/******************************************
混合正規分布を構成する多変量正規分布の
パラメータを格納する構造体
******************************************/
typedef struct {
size_t ndim; // 次元数
double *mu; // 平均ベクトル（ndim 次元）
double *cov; // 分散共分散行列（ndim(行) * ndim(列）次元)
double detcov; // 分散共分散行列の行列式
double *icov; // 分散共分散行列の逆行列
} gaussian_params;
/*** 正規分布パラメータ構造体の領域確保 ***/
gaussian_params *
gaussian_params_alloc (void)
{
gaussian_params *par = (gaussian_params *) malloc (sizeof (gaussian_params));
par->ndim = -1;
par->mu = NULL;
par->cov = NULL;
par->icov = NULL;
return par;
}
/*** 正規分布パラメータ構造体の開放 ***/
void
gaussian_params_free (gaussian_params *par)
{
if (par) {
if (par->mu) free (par->mu);
if (par->cov) free (par->cov);
if (par->icov) free (par->icov);
free (par);
}
return;
}
構造体の名前は gaussian_params としました。
で、この中の行列式や逆行列は GSL の関数群を使って計算する事にしましょう。
ただ、ユーザがそれらを計算して各々代入するのも手間なので、共分散行列を代入済みの gaussian_params 構造体を引数に、これらを勝手に計算させる関数を用意します。
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
/****************************************************
LU分解を用い共分散行列の行列式・逆行列を計算
gaussian_params のメンバ double *cov から
行列式と逆行列を計算し、メンバ double detcov
と double *icov に計算結果をそれぞれ格納する。
****************************************************/
void
gaussian_params_det_and_inv_covariance (gaussian_params *par)
{
int s;
// Permutation
gsl_permutation *p;
gsl_matrix *lu;
// 1次元配列 par->cov に対する行列の像
// Matrix_view
gsl_matrix_view cov;
// 1次元配列 par->icov に対する行列の像
gsl_matrix_view icov;
cov = gsl_matrix_view_array (par->cov, par->ndim, par->ndim);
if (par->icov == NULL)
par->icov = (double *) malloc (par->ndim * par->ndim * sizeof (double));
icov = gsl_matrix_view_array (par->icov, par->ndim, par->ndim);
p = gsl_permutation_alloc (par->ndim);
// LU分解を計算するために用いるテンポラリな行列 *lu に cov.matrix をコピー
lu = gsl_matrix_alloc (par->ndim, par->ndim);
gsl_matrix_memcpy (lu, &cov.matrix);
// LU-Decomposition
gsl_linalg_LU_decomp (lu, p, &s); // LU分解
par->detcov = gsl_linalg_LU_det (lu, s); // 行列式
gsl_linalg_LU_invert (lu, p, &icov.matrix); // 逆行列
gsl_permutation_free (p);
gsl_matrix_free (lu);
return;
}
また、構造体 gaussian_params を作成する関数として以下を用意します。
ユーザが指定した次元数、平均、共分散行列をコピーして保持します。
ポインタでも良い様な気もしますが、EMの更新で書き換えられるのでコピーして初期値とするようにしました。
またこの中で行列式・逆行列も計算します。
/********************************************************
gaussian_params 構造体を作成してそのポインタを返す
引数 :
size_t ndim : 次元数
double *mu : 平均ベクトルの初期値
double *var : 共分散行列の初期値
引数で与えられた上記を、それぞれメンバ *mu, *var に
コピーする。
********************************************************/
gaussian_params *
gaussian_params_new (size_t ndim, double *mu, double *cov)
{
int j;
gaussian_params *par;
par = gaussian_params_alloc ();
// 次元数を格納
par->ndim = ndim;
// 入力された平均ベクトルをコピー
par->mu = (double *) malloc (ndim * sizeof (double));
for (j = 0; j < ndim; j++) par->mu[j] = mu[j];
// 入力された共分散行列をコピー
par->cov = (double *) malloc (ndim * ndim * sizeof (double));
for (j = 0; j < ndim * ndim; j++) par->cov[j] = cov[j];
// 逆行列、行列式を計算
par->icov = (double *) malloc (ndim * ndim * sizeof (double));
gaussian_params_det_and_inv_covariance (par);
return par;
}
これらの関数は機能が少し多すぎな気もしますが、あんまりいっぱい関数をつくるのもブログに載せるには見づらいような気がしますので。。。
あとエラー処理が（メンドクサイノデ）なんにも書かれてませんが（汗
おいおい追加していきましょう（絶対しないな・・・）
まぁ、こんなところをきっかけとして、次回以降EMアルゴリズムのプログラムを作って逝きましょう。