ifa.cpp

变分独立因子分析C++代码,H.Attias说比独立分量分析要好,但是这个程序分析效果不好,可能是程序问题,也可能是对理论理解不透

// W.cpp: implementation of the W class.
//
//////////////////////////////////////////////////////////////////////

#include "ifa.h"
#include <cmath>
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

const double pi2=6.283185307179586200;

W::W()
{
	Lp=0;
	L=0;
	maxni=0;
	totalni=0;
	int* n=0;
}

W::~W()
{

}


void W::initialize(int l, int lp, int *ni)
{
	L=l;
	Lp=lp;
	n=new int[L];
	totalni=0;
	maxni=0;
	for(int i=0;i<L;i++) 
	{
		n[i]=ni[i];
		totalni+=n[i];
		if(maxni<n[i]) maxni=n[i];
	}

	H=CMatrix(Lp,L);
	for(i=0;i<Lp;i++)
	{
		for(int j=0;j<L;j++)
		{
			H[i][j]=1.0/L;
		}
	}

	lambda=CMatrix(Lp,Lp);
	for(i=0;i<Lp;i++)
	{
		for(int j=0;j<Lp;j++)
		{
			lambda[i][j]=0;
		}
		lambda[i][i]=1;
	}

	mu=CMatrix(L, maxni);
	for(i=0;i<L;i++)
	{
		for(int j=0;j<n[i];j++) 
		{
			mu[i][j]=1.0/n[i];
		}
	}

	nu=CMatrix(L,maxni);
	for(i=0;i<L;i++)
	{
		for(int j=0;j<n[i];j++)
		{
			nu[i][j]=0;
		}
	}

	omega=CMatrix(L,maxni);
	for(i=0;i<L;i++)
	{
		for(int j=0;j<n[i];j++)
		{
			omega[i][j]=1.0/n[i];
		}
	}
}

double pq(const W & w, const valarray<int> & q)
{
	double pq=1;
	for(int i=0;i<w.L;i++)
	{
		pq*=w.omega[i][q[i]];
	}
	return pq;
}

double py_q(const W & w, const CMatrix & y, const valarray<int> & q)
{
	CMatrix muq(w.L,1);
	for(int i=0;i<w.L;i++)
	{
		muq[i][0]=w.mu[i][q[i]];
	}

	CMatrix Vq(w.L,w.L);
	for(i=0;i<w.L;i++)
	{
		for(int j=0;j<w.L;j++)
		{
			Vq[i][j]=0;
		}
		Vq[i][i]=w.nu[i][q[i]];
	}

	CMatrix H=w.H;
	CMatrix lambda=w.lambda;
	return Gauss(y, H*muq,H*Vq*(~H)+lambda);
}

double py(const W & w, const CMatrix & y)
{
	double py=0;
	valarray<int>  q(0,w.L);
	int i=0;
	while(i>=0)
	{
		if(i==w.L)
		{
			py+=pq(w,q)*py_q(w,y,q);
			i--;
			q[i]++;
		}
		else
		{
			if(q[i]<w.n[i])
			{
				i++;
			}
			else
			{
				q[i]=0;
				i--;// if(i<0) break;
				if(i>=0) q[i]++;// when i=-1，the sentence q[i]++ will make i be 0
			}
		}		
	}
	return py;	
}

double pq_y(const W & w, const valarray<int> & q, const CMatrix & y, const double py)
{
	return pq(w, q)*py_q(w,y,q)/py;
}

double Gauss(const CMatrix & x, const CMatrix & mu, const CMatrix & sigmma)
{
	CMatrix xp=x;
	CMatrix mup=mu;
	CMatrix sigmmap=sigmma;
	return 1/sqrt(symmdeterminant(sigmmap*pi2))\
		*exp(-0.5*scalarvalue( (~(xp-mup))*symminverse(sigmmap)*(xp-mup) ));
}

CMatrix sigmaq(const W & w, const valarray<int> & q)
{
	CMatrix Vq(w.L,w.L);
	for(int i=0;i<w.L;i++)
	{
		for(int j=0;j<w.L;j++)
		{
			Vq[i][j]=0;
		}
		Vq[i][i]=w.nu[i][q[i]];
	}

	return symminverse( ~w.H*symminverse(w.lambda)*w.H + diaginverse(Vq) );

}

CMatrix rhoqy(const W & w, const valarray<int> & q, const CMatrix & y,  CMatrix & sigmaq)//const
{
	CMatrix Vq(w.L,w.L);
	for(int i=0;i<w.L;i++)
	{
		for(int j=0;j<w.L;j++)
		{
			Vq[i][j]=0;
		}
		Vq[i][i]=w.nu[i][q[i]];
	}
	
	CMatrix muq(w.L,1);
	for( i=0;i<w.L;i++)
	{
		muq[i][0]=w.mu[i][q[i]];
	}

	//CMatrix sigma=sigmaq;

	return sigmaq*( (~w.H)*symminverse(w.lambda)*y + diaginverse(Vq)*muq );
}


CMatrix Attias(const valarray<int> & N, const CMatrix & Vq, const CMatrix &H, \
			   const CMatrix & lambda, const CMatrix & y, const CMatrix & muq)
{
	CMatrix psi(N.size(),N.max());
	return psi;
}