まずはユーティリティ群を。
#include <gsl gsl_math.h>
#include <gsl gsl_rng.h>
#include <gsl gsl_randist.h>
const gsl_rng *_rng_;
/*** ユーティリティー ***/
/* ベクトルの和をとる */
double
_sum (int n, double *vector)
{
int i;
double sum = 0.;
for (i = 0; i < n; i++) sum += vector[i];
return sum;
}
/* ベクトルをその和で規格化する */
void
_normalize (int n, double *vector)
{
int i;
double sum = _sum (n, vector);
for (i = 0; i < n; i++) vector[i] /= sum;
return;
}
/* ベクトルを出力 */
void
fprintf_vector (FILE *stream, int n, double *vector, char *format)
{
int i;
for (i = 0; i < n; i++) {
fprintf (stream, format, vector[i]);
if (i < n - 1) fprintf (stream, "\t");
}
fprintf (stream, "\n");
return;
}
/* 行列を出力 */
void
fprintf_matrix (FILE *stream, int nrow, int ncol, double **matrix, char *format)
{
int i;
for (i = 0; i < nrow; i++) fprintf_vector (stream, ncol, matrix[i], format);
return;
}
これらでベクトル(一次元配列)の和の計算、ベクトルをその和で規格化、ベクトルの出力及び行列（2次元配列）の出力を行います。
次にEMアルゴリズムの部分。
/*** EMアルゴリズム ***/
/* 対数尤度の条件付き期待値 */
double
log_likelihood (int ncls, int nitems, double *xi, double **z, int **x, double **theta)
{
int i, k;
double ll = 0.;
// ll = sum_i sum_k { z[i][k] * log (xi[k] * f (x[i] | theta[k]) ) }
for (i = 0; i < nitems; i++) {
for (k = 0; k < ncls; k++) {
ll += z[i][k] * (gsl_ran_multinomial_lnpdf (nitems, theta[k], x[i]) + log (xi[k]));
}
}
return ll;
}
/* Eステップ */
void
E_step (int n, int ncls, int nitems, double *xi, double **z, int **x, double **theta)
{
int i, k;
for (i = 0; i < n; i++) {
for (k = 0; k < ncls; k++) {
// z[i][k] = xi[k] * f(x[i] | theta[k])
z[i][k] = xi[k] * gsl_ran_multinomial_pdf (nitems, theta[k], x[i]);
}
// z[i][k] = z[i][k] / sum_k z[i][k]
_normalize (ncls, z[i]);
}
return;
}
/* Mステップ */
void
M_step (int n, int ncls, int nitems, double *xi, double **z, int **x, double **theta)
{
int i, k, l;
for (k = 0; k < ncls; k++) {
/* theta の更新 */
// theta[k][l] = sum_i (z[i][k] * x[i][l])
for (l = 0; l < nitems; l++) {
theta[k][l] = 0.;
for (i = 0; i < n; i++) theta[k][l] += z[i][k] * ((double) x[i][l]);
}
// theta[k] = theta[k] / sum_l theta[k][l]
_normalize (nitems, theta[k]);
/* xi の更新 */
// xi[k] = sum_i z[i][k] / n
xi[k] = 0.;
for (i = 0; i < n; i++) xi[k] += z[i][k];
xi[k] /= (double) n;
}
return;
}
一番目の関数で対数尤度の条件付期待値を計算します。この時、gsl の関数
double gsl_ran_multinomial_lnpdf (size_t K, const double p[], const unsigned int x[]);
という関数を使っています。
これはパラメータ const double p[] （次元は size_t K）の多項分布でデータ const unsigned int x[] が得られる密度関数の log を計算します。
ちなみに
double gsl_ran_multinomial_pdf (size_t K, const double p[], const unsigned int x[]);
が密度関数、
void gsl_ran_multinomial_lnpdf (const gsl_rng *r, size_t K, unsigned int n, const double p[], unsigned int x[]);
が多項分布の乱数発生器になります。
つづいて入力とするデータを作成する部分
/*** テストデータ作成用 ***/
/* パラメータを作成 (nitems 次元ベクトルを ncls 個分) */
double **
create_parameter (int ncls, int nitems)
{
int k, l;
double **theta = (double **) malloc (ncls * sizeof (double *));
for (k = 0; k < ncls; k++) {
theta[k] = (double *) malloc (nitems * sizeof (double));
for (l = 0; l < nitems; l++) theta[k][l] = gsl_ran_flat (_rng_, 0., 1.);
_normalize (nitems, theta[k]);
}
return theta;
}
/* データ作成 (nitems 次元ベクトルを n 個分) */
int **
create_data (int n, int ncls, int nitems, int *m, double *xi, double **theta)
{
int i, k;
int **x;
x = (int **) malloc (n * sizeof (int *));
for (i = 0; i < n; i++) {
// 混合率 xi に応じて多項分布(パラメータを指定する k の値)を選ぶ
double r = gsl_ran_flat (_rng_, 0., 1.);
for (k = 0; k < ncls; k++) {
if (r < xi[k]) break;
r -= xi[k];
}
x[i] = (int *) malloc (nitems * sizeof (int));
gsl_ran_multinomial (_rng_, nitems, m[i], theta[k], x[i]);
}
return x;
}
その際、各々の多項分布のパラメータの和が 1 になるよう規格化しています。
2番目の関数で混合多項分布に従うデータを作成しています。
与えられた混合比に応じてパラメータを選び、乱数を生成して2次元配列に入れて返します。
最後にメイン関数。
int
main (void)
{
int ncls = 3; // コンポーネント数 (K)
int nitems = 4; // アイテム数 (L)
int n = 1000; // データセット数 (N)
int **x; // データ
double **theta; // パラメータ
double *xi; // 混合比
double **z; // 事後確率
_rng_ = gsl_rng_alloc (gsl_rng_default); // gslの乱数初期化
// 混合多項分布に従うデータの作成
{
int i;
double **theta0 = create_parameter (ncls, nitems); // 各多項分布のパラメータ
double xi0[3] = {0.2, 0.5, 0.3}; // 混合率
int *m0; // 試行回数 ( m0[i] = sum_l x[i][l] )
m0 = (int *) malloc (n * sizeof (int));
for (i = 0; i < n; i++) m0[i] = 500 + gsl_rng_uniform_int (_rng_, 100); // 500～600回
x = create_data (n, ncls, nitems, m0, xi0, theta0);
fprintf (stdout, "xi0 :\n");
fprintf_vector (stdout, ncls, xi0, "%.3f");
fprintf (stdout, "\ntheta0 :\n");
fprintf_matrix (stdout, ncls, nitems, theta0, "%.3f");
}
// 領域確保と初期値設定
{
int i, k, l;
theta = (double **) malloc (ncls * sizeof (double *));
for (k = 0; k < ncls; k++) {
theta[k] = (double *) malloc (nitems * sizeof (double));
for (l = 0; l < nitems; l++) theta[k][l] = 1. / (double) nitems;
}
z = (double **) malloc (n * sizeof (double *));
for (i = 0; i < n; i++) {
z[i] = (double *) malloc (ncls * sizeof (double));
for (k = 0; k < ncls; k++) z[i][k] = 0.;
}
xi = (double *) malloc (ncls * sizeof (double));
for (k = 0; k < ncls; k++) xi[k] = gsl_ran_flat (_rng_, 0., 1.);
_normalize (ncls, xi);
}
// EMアルゴリズムを実行
{
int iter = 0;
double l_prev;
double l = -GSL_DBL_MAX;
while (1) {
l_prev = l;
E_step (n, ncls, nitems, xi, z, x, theta);
l = log_likelihood (ncls, nitems, xi, z, x, theta);
if (fabs (l - l_prev) < 1.e-16 || ++iter >= 1000) break;
M_step (n, ncls, nitems, xi, z, x, theta);
}
fprintf (stdout, "\niter = %d, l = %f\n\n", iter, l);
fprintf (stdout, "xi :\n");
fprintf_vector (stdout, ncls, xi, "%.3f");
fprintf (stdout, "\ntheta :\n");
fprintf_matrix (stdout, ncls, nitems, theta, "%.3f");
}
return 0;
}
<真のパラメータ>
xi0 :
0.200 0.500 0.300
theta0 :
0.418 0.068 0.118 0.396
0.096 0.201 0.396 0.308
0.200 0.274 0.282 0.244
iter = 15, l = -41.475837
<計算結果>
xi :
0.498 0.317 0.185
theta :
0.096 0.200 0.396 0.309
0.201 0.275 0.281 0.243
0.419 0.066 0.119 0.396
となりました。
でも少し条件を変えると収束が遅くなったり、更新されなくなったりします。
そういうものなのか、或いは何処か間違ってるのか・・・多分後者だと思う（汗