やっと基本のEMアルゴリズム

これまで勾配法で最適解を求めてきましたが、今回はEMアルゴリズムを使って解を求めていきます。

EMアルゴリズムとは、E（Expectation）ステップとM（Maximization）ステップを繰り返していき、解を求めるアルゴリズムです。

今回は、ガウス混合分布をですとデータとしてパラメータを推定します。

まずデータとして、
$N(x_{1}|1,1),N(x_{2}|-1,1)$ の正規分布にのっとったデータを生成します。この時、 $p(x_{1})=0.4,p(x_{2})=0.6$ としました。

この時の混合ガウス分布は

$p({\bf x})=\Sigma^{K}_{k=1}\pi_{k}N({\bf x}|{\bf \mu_{k},\Sigma_{k}})$

となります。

今回求めたいのは、 $\pi_{k}$ です。ただし、
$0 \leq \pi_{k} \leq 1 , \Sigma^{K}_{k=1}\pi_{k} = 1$
という条件付きです。

基本的な流れは

１．平均・分散（今回は１に設定して省略）・混合係数 $\pi_{k}$ を適当に設定

２．現在のパラメータから事後確率計算
＊事後確率はデータ数×混合係数の数(パラメータの数？)だけ発生する。

３．事後確率を用いてパラメータ更新

４．対数尤度を計算し、変化が小さくなれば終了。対数尤度は

$ln({\bf X}|{\bf \mu,\Sigma, \pi})=\Sigma^{N}_{n=1}ln\{\Sigma^{K}_{k=1}\pi_{k}N({\bf x}|{\bf \mu}_{k},{\bf \Sigma}_{k})\}$

となります。

本当は事後確率のあたりの式をしっかり描くべきですが、もう少し読み込まないと自信がないので、今回はこんな感じで行きます。（理論ありきではあるのですが、コーディングから入るのも悪くないかなと自分では思います。）

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <vector>
#include <algorithm>

using namespace std;
int Dm,Dn,Dk;
double **MatAverage,**MatData,**MatGamma,*VecTheta;

void readValue(char *fn1){
  ifstream fin;
  int i,j,k,l;
  double v;
  fin.open(fn1);
  fin >> Dm >> Dn >> k;
  MatData = new double*[Dm];
  for(i=0;i<Dm;i++){
    MatData[i] = new double[Dn];
    for(j=0;j<Dn;j++){
      fin >> MatData[i][j];
    }
    fin >> k >> l;
  } 
}

void initValue(){
  int i,j,k;
  MatAverage = new double *[Dk];
  MatGamma = new double *[Dm];
  VecTheta = new double [Dk];
  for(k=0;k<Dk;k++) VecTheta[k] = 1.0/Dk;
  for(i=0;i<Dm;i++) MatGamma[i] = new double[Dk];
  for(k=0;k<Dk;k++){
    MatAverage[k] = new double[Dn];
    for(j=0;j<Dn;j++) MatAverage[k][j] = ((rand()%20000)/10000.0)-1.0;
  }
}

double Expectation(){
  int i, j, k;
  double v, w, x;
  for(i=0,w=0.0;i<Dm;i++){
    for(k=0,v=0.0;k<Dk;k++){
      for(j=0,x=0.0;j<Dn;j++){
	x += (MatData[i][j] - MatAverage[k][j])*(MatData[i][j] - MatAverage[k][j]);
      }
      MatGamma[i][k] = VecTheta[k]*exp(-0.5*x);
      v += MatGamma[i][k];
    }
    for(k = 0; k < Dk; k++) MatGamma[i][k] /= v;
    w += log(v);
  }
  return(w);
}

void Maximization(){
  int i, j, k;
  for(k=0;k<Dk;k++){
    VecTheta[k] = 0.0;
    for(j=0;j<Dn;j++) MatAverage[k][j] = 0.0;
  }
  for(i=0;i<Dm;i++){
    for(k=0;k<Dk;k++){
      VecTheta[k] += MatGamma[i][k];
      for(j=0;j<Dn;j++) MatAverage[k][j] += MatGamma[i][k] * MatData[i][j];
    }
  }
  for(k=0;k<Dk;k++){
    for(j=0;j<Dn;j++) MatAverage[k][j] /= VecTheta[k];
    VecTheta[k] /= Dm;
  }
}

void printValue(char *fn1){
  int j, k;
  for(k = 0; k < Dk; k++) printf("%e ", VecTheta[k]);
  printf("\n");
  for(k = 0; k < Dk; k++){
    for(j = 0; j < Dn; j++) printf("%e ", MatAverage[k][j]);
    printf("\n");
  }
}

int main(int argc, char **argv){
  int i,j,k;
  double Likelihood,oldL=0.0;
  srand(time(0));
  Dk = atoi(argv[1]); //set param
  readValue(argv[2]);
  initValue();
  printf("init ok\n");
  for(i=0;i<1000;i++){
    Likelihood = Expectation();
    printf("%d %e\n",i+1,Likelihood);
    Maximization();
    if(fabs(Likelihood - oldL) < 1.0e-8) break;
    oldL = Likelihood;
  }
  printValue(argv[3]);
  printf("print ok\n");
}

結果

媒介変数
3.926149e-01 6.073851e-01

平均
1.070515e+00 -9.543940e-01

さらに、対数尤度のプロットです。このプロットはid:syou6162さんの
[http://d.hatena.ne.jp/syou6162/20091018/1255826530:[R][機械学習]初めてのEMアルゴリズム with R]
を参考にさせていただきました。

尤度関数の増加がわかります。

まだまだ理論的な部分の理解が必要です。また、これが理解できないと今の研究で生きていけないので、もっと精進したいと思います。