em.cpp

用C++实现了期望最大化算法

// Em.cpp: implementation of the CEm class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "imageassociation.h"
#include "Em.h"
#include <math.h>
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
template<class T>
CEm<T>::CEm()
{

}
template<class T>
CEm<T>::~CEm()
{

}
template<class T>
void CEm<T>::initialize(T **data,int num,int dim)//double
{
  m_Data  = data;
  m_N = num;//Data.length;
  m_d = dim;//Data[0].length;
  m_Min   = new double[m_d];
  m_Max   = new double[m_d];
  int n,i;//,j;
  for(i=0;i<m_d;i++) {
    m_Min[i] = m_Data[0][i];
    m_Max[i] = m_Data[0][i];
    for(n=1;n<m_N;n++) {
      if(m_Min[i] > m_Data[n][i]) m_Min[i] = m_Data[n][i]; else
      if(m_Max[i] < m_Data[n][i]) m_Max[i] = m_Data[n][i];
    }
  }
}
template<class T>
void CEm<T>::setParameter(int m)
{
  int i,j;
  double var = 0.0;
  m_M = m;
  m_P = new double[m];
  m_u = new  double*[m];
  for(i=0;i<m;i++){// [m_d];
	  m_u[i]=new double[m_d];
  }
  m_v = new double[m];
  m_P_old = new double[m];
  m_u_old=new double*[m];//
  for(i=0;i<m;i++){
	  m_u_old[m]=new double[m_d];
  }
  m_v_old = new double[m];
  m_P_x   = new double[m];
  m_P_x_old = new double[m_N];
  for(i=0;i<m_d;i++)
    var += (m_Max[i] - m_Min[i]) * (m_Max[i] - m_Min[i]) / 12 / m_d;
  for(j=0;j<m;j++)
  {
    m_P[j] = 1.0 / (double)m;
    m_v[j] = var / 4 + random() * var / 4; 
    for(i=0;i<m_d;i++)
      m_u[j][i] = (m_Max[i] + m_Min[i]) / 2 +
                (random() - 0.5 ) *  (m_Max[i] - m_Min[i]) / 2;
  }
}
template<class T>
void CEm<T>::new_old()
{
  int i,j;
  for(j=0;j<m_M;j++)
  {
    m_P_old[j] = m_P[j];
    m_v_old[j] = m_v[j];
    for(i=0;i<m_d;i++)
		m_u_old[j][i] = m_u[j][i];
  }
}
template<class T>
double CEm<T>::P_x_j(int n, int j)
{
  int i;
  double e = 0.0;
  for(i=0;i<m_d;i++)
    e += (m_Data[n][i] - m_u[j][i]) * (m_Data[n][i] - m_u[j][i]);
  e /= -2 * m_v[j] * m_v[j];
  return exp(e) / pow(2 * PI * m_v[j] * m_v[j], (double)m_d / 2.0); 
} 

 /* double CEm::P_j_x(int n, int j)
{
  return P_x_j(n,j) * P[j] / P_x[n];
}//*/
template<class T>
void CEm<T>::P_x()
{
  int n,j;
  for(n=0;n<m_N;n++)
  {
    double temp = 0.0;
    for(j=0;j<m_M;j++)
      temp += P_x_j(n,j) * m_P[j];
    m_P_x[n] = temp;
  }
}
template<class T>
double CEm<T>::P_x_j_old(int n, int j)
{
  int i;
  double e = 0.0;
  for(i=0;i<m_d;i++)
    e += (m_Data[n][i] - m_u_old[j][i]) * (m_Data[n][i] - m_u_old[j][i]);
  e /= -2 * m_v_old[j] * m_v_old[j];
  return exp(e) / pow(2 * PI * m_v_old[j] * m_v_old[j], (double)m_d / 2.0); 
} 
template<class T>
double CEm<T>::pdf(double x[])
{
  int i,j;
  double e = 0.0;
  double t = 0.0;
  for(j=0;j<m_M;j++) {
    for(i=0;i<m_d;i++)
      e += (x[i] - m_u[j][i]) * (x[i] - m_u[j][i]);
    e /= -2 * m_v[j] * m_v[j];
    e = exp(e) / pow(2 * PI * m_v[j] * m_v[j],(double)m_d / 2.0); 
    t += e*m_P[j];
  }
  return t;
}

/*double CEm::P_j_x_old(int n, int j)
{
  return P_x_j_old(n,j) * P_old[j] / P_x_old[n];
}*/
template<class T>
void CEm<T>::P_x_old()
{
  int n,j;
  for(n=0;n<m_N;n++)
  {
    double temp = 0.0;
    for(j=0;j<m_M;j++)
      temp += P_x_j_old(n,j) * m_P_old[j];
    m_P_x_old[n] = temp;
  }
}
template<class T>
void CEm<T>::iteration()
{
	  int i,j,n;
	  new_old();
	  P_x();
	  P_x_old();
	  
	  // calculate new mean & variance
	  
	  for(j=0;j<m_M;j++)
	  {
		  double denom = 0.0;
		  for(n=0;n<m_N;n++)
			  denom += P_j_x_old(n,j);
		  
		  for(i=0;i<m_d;i++)
		  { 
			  m_u[j][i] = 0.0;
			  for(n=0;n<m_N;n++)
				  m_u[j][i] += P_j_x_old(n,j) * m_Data[n][i];
			  m_u[j][i] /= denom;
		  }
		  
		  m_v[j] = 0.0;
		  for(n=0;n<m_N;n++)
		  {
			  double sum = 0.0;
			  for(i=0;i<m_d;i++)
				  sum += (m_Data[n][i] - m_u[j][i]) * (m_Data[n][i] - m_u[j][i]);
			  m_v[j] += sum * P_j_x_old(n,j);
		  }
		  m_v[j] /= denom * m_d;
		  m_v[j] = sqrt(m_v[j]);
		  
		  m_P[j] = denom / (double)m_N;
	  }
}


template<class T>
void CEm<T>::display(HDC pDC)
{
	  CString newline;// = System.getProperty("line.separator");
	  int i,j;
	  for(j=0;j<m_M;j++)
	  {
		  newline.Format("P(%d)=%f	V(%d)=%f	U(%d)",j,m_P[j]);
		  CString temp;
		  for(i=0;i<m_d;i++){
			  //System.out.print(u[j][i] + " ");
			  temp.Format("u(%d)(%d)=%f",j,i,m_u[j][i]);
			  newline+=temp;
		  }
		  TextOut(pDC,0,j*15,newline,newline.GetLength());
		  
	  }
}