prepare structure for R package

[valse.git] / OLD_MATLAB / ProcLassoRank / EMGrank.m

function[phi,LLF] = EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank)\r
\r
	% get matrix sizes\r
	[~,m,k] = size(Rho);\r
	[n,p] = size(X);\r
\r
	% allocate output matrices\r
	phi = zeros(p,m,k);\r
	Z = ones(n,1,'int64');\r
	LLF = 0.0;\r
\r
	% local variables\r
	Phi = zeros(p,m,k);\r
	deltaPhi = 0.0;\r
	deltaPhi = [];\r
	sumDeltaPhi = 0.0;\r
	deltaPhiBufferSize = 20;\r
\r
	%main loop (at least mini iterations)\r
	ite = int64(1);\r
	while ite<=mini || (ite<=maxi && sumDeltaPhi>tau)\r
\r
		%M step: Mise à jour de Beta (et donc phi)\r
		for r=1:k\r
			if (sum(Z==r) == 0)\r
				continue;\r
			end\r
			%U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr\r
			[U,S,V] = svd(pinv(transpose(X(Z==r,:))*X(Z==r,:))*transpose(X(Z==r,:))*Y(Z==r,:));\r
			%Set m-rank(r) singular values to zero, and recompose \r
			%best rank(r) approximation of the initial product\r
			S(rank(r)+1:end,:) = 0;\r
			phi(:,:,r) = U * S * transpose(V) * Rho(:,:,r);\r
		end\r
		\r
		%Etape E et calcul de LLF\r
		sumLogLLF2 = 0.0;\r
		for i=1:n\r
			sumLLF1 = 0.0;\r
			maxLogGamIR = -Inf;\r
			for r=1:k\r
				dotProduct = (Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r));\r
				logGamIR = log(Pi(r)) + log(det(Rho(:,:,r))) - 0.5*dotProduct;\r
				%Z(i) = index of max (gam(i,:))\r
				if logGamIR > maxLogGamIR\r
					Z(i) = r;\r
					maxLogGamIR = logGamIR;\r
				end\r
				sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2);\r
			end\r
			sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
		end\r
		\r
		LLF = -1/n * sumLogLLF2;\r
\r
		% update distance parameter to check algorithm convergence (delta(phi, Phi))\r
		deltaPhi = [ deltaPhi, max(max(max((abs(phi-Phi))./(1+abs(phi))))) ];\r
		if length(deltaPhi) > deltaPhiBufferSize\r
			deltaPhi = deltaPhi(2:length(deltaPhi));\r
		end\r
		sumDeltaPhi = sum(abs(deltaPhi));\r
\r
		% update other local variables\r
		Phi = phi;\r
		ite = ite+1;\r
\r
	end\r
\r
end\r