em.h

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

// Em.h: interface for the CEm class.
//
//////////////////////////////////////////////////////////////////////
// 假设了各维相互独立，只有方差，没有协方差
#if !defined _EM_H_
#define _EM_H_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#ifndef PI
#define PI 3.141592653589793115997963468544185161590576171875
#endif
template<class T>//,int num,int dim>
class CEm  
{
public:
	T **	   m_Data;//double
	int        m_M;
	int        m_N;
	int        m_d;
	
	double *   m_P;		// Mixing Parameters
	double **  m_u;		// Mean
	double *   m_v;     // Variance
	
	double *   m_P_old;	// Mixing Parameters
	double **  m_u_old;	// Mean
	double *   m_v_old; // Variance
	
	double *   m_P_x;	// P_x calculation results
	double *   m_P_x_old;	// P_x calculation results
	
	T *        m_Min;//double	// Maximum of each dimension
	T *        m_Max;//double	// Minimum of each dimension
	CEm();
	virtual ~CEm();
	CEm(T **data,int num,int dim, int m) {//double
		initialize(data,num,dim);
		setParameter(m);
	};
	
	CEm(T **data,int num,int dim) {//double
		initialize(data,num,dim);
	};
    void initialize(T **data,int num,int dim);//double
	void setParameter(int m);
	void new_old();
	double P_x_j(int n, int j);
	double P_j_x(int n, int j)  { return P_x_j(n,j) * m_P[j] / m_P_x[n]; };
	double P_x_j_old(int n, int j);
	double pdf(double x[]);
	double P_j_x_old(int n, int j){ return P_x_j_old(n,j) * m_P_old[j] / m_P_x_old[n]; } ;
	void P_x();
	void P_x_old();
	void iteration();
	void display(HDC pDC);
};
inline double random(){
	//srand( (unsigned)time( NULL ) ); //有这句反而产生相同的数
	double result=(double)rand()/RAND_MAX;
	ASSERT(result>=0&&result<=1);
	return result;
}
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 T[m_d];
  m_Max   = new T[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[i]=new double[m_d];
  }
  m_v_old = new double[m];
  m_P_x   = new double[m_N];
  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(int n=0;n<m_N;++n)//
//	  for(i=0;i<m_d;++i)//
//		  var+=m_Data[n][i]*m_Data[n][i];//
//  var/=m_N;//
  for(j=0;j<m;j++)
  {
    m_P[j] = 1.0 / (double)m;
    m_v[j] = var / 4 + random() * var / 4; 
	//m_v[j]=var;
	//n=j*(m_N-1)/(m-1);
	//for(i=0;i<m_d;i++)
	//	m_u[j][i]=m_Data[n][i];
    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>::setParameter(int m, T* means, T** delta)// 假设了各维相互独立，只有方差，没有协方差
{
	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_u[i]=new double[m_d];
	}
	for(i=0;i<m;++i){
		for(j=0;j<m_d;++j){
			m_u[i][j]=means[i*m_d+j];
		}
	}
	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[i]=new double[m_d];
	}
	m_v_old = new double[m];
	m_P_x   = new double[m_N];
	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(int n=0;n<m_N;++n)//
	//	  for(i=0;i<m_d;++i)//
	//		  var+=m_Data[n][i]*m_Data[n][i];//
	//  var/=m_N;//
	for(j=0;j<m;j++)
	{
		m_P[j] = 1.0 / (double)m;
		m_v[j] = var / 4 + random() * var / 4; 
		//m_v[j]=var;
		//n=j*(m_N-1)/(m-1);
		//for(i=0;i<m_d;i++)
		//	m_u[j][i]=m_Data[n][i];
		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], 1 / 2);
	  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], 1 / 2); 
	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 ",j,m_P[j],j,m_v[j]);
		  CString temp;
		  for(i=0;i<m_d;i++){
			  temp.Format("u(%d)(%d)=%f ",j,i,m_u[j][i]);
			  newline+=temp;
		  }
		  TextOut(pDC,0,j*15,newline,newline.GetLength());
	  }
}
#endif // !defined _EM_H_

/*
    //CEm使用实例
	CEm<float> em;
	em.initialize(data,num,dim);
	em.setParameter(5);//models=10.
	for(i=0;i<10;++i){
		//em.new_old();
		em.iteration();
	}
	em.display(pDC->m_hDC);
*/