initial commit

[valse.git] / ProcLassoMLE / EMGLLF.m

function[phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)\r
\r
	%Get matrices dimensions\r
	PI = 4.0 * atan(1.0);\r
	n = size(X, 1);\r
	[p,m,k] = size(phiInit);\r
	\r
	%Initialize outputs\r
	phi = phiInit;\r
	rho = rhoInit;\r
	pi = piInit;\r
	LLF = zeros(maxi,1);\r
	S = zeros(p,m,k);\r
	\r
	%Other local variables\r
	%NOTE: variables order is always n,p,m,k\r
	gam = gamInit;\r
	Gram2 = zeros(p,p,k);\r
	ps2 = zeros(p,m,k);\r
	b = zeros(k,1);\r
	pen = zeros(maxi,k);\r
	X2 = zeros(n,p,k);\r
	Y2 = zeros(n,m,k);\r
	dist = 0;\r
	dist2 = 0;\r
	ite = 1;\r
	pi2 = zeros(k,1);\r
	ps = zeros(m,k);\r
	nY2 = zeros(m,k);\r
	ps1 = zeros(n,m,k);\r
	nY21 = zeros(n,m,k);\r
	Gam = zeros(n,k);\r
	EPS = 1e-15;\r
\r
	while ite<=mini || (ite<=maxi && (dist>=tau || dist2>=sqrt(tau)))\r
		\r
		Phi = phi;\r
		Rho = rho;\r
		Pi = pi;\r
		\r
		%Calculs associés à Y et X\r
		for r=1:k\r
			for mm=1:m\r
				Y2(:,mm,r) = sqrt(gam(:,r)) .* Y(:,mm);\r
			end\r
			for i=1:n\r
				X2(i,:,r) = X(i,:) .* sqrt(gam(i,r));\r
			end\r
			for mm=1:m\r
				ps2(:,mm,r) = transpose(X2(:,:,r)) * Y2(:,mm,r);\r
			end\r
			for j=1:p\r
				for s=1:p\r
					Gram2(j,s,r) = dot(X2(:,j,r), X2(:,s,r));\r
				end\r
			end\r
		end\r
\r
		%%%%%%%%%%\r
		%Etape M %\r
		%%%%%%%%%%\r
		\r
		%Pour pi\r
		for r=1:k\r
			b(r) = sum(sum(abs(phi(:,:,r))));\r
		end\r
		gam2 = sum(gam,1);\r
		a = sum(gam*transpose(log(pi)));\r
		\r
		%tant que les proportions sont negatives\r
		kk = 0;\r
		pi2AllPositive = false;\r
		while ~pi2AllPositive\r
			pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi);\r
			pi2AllPositive = true;\r
			for r=1:k\r
				if pi2(r) < 0\r
					pi2AllPositive = false;\r
					break;\r
				end\r
			end\r
			kk = kk+1;\r
		end\r
		\r
		%t(m) la plus grande valeur dans la grille O.1^k tel que ce soit\r
		%décroissante ou constante\r
		while (-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*transpose(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000\r
			pi2 = pi+0.1^kk*(1/n*gam2-pi);\r
			kk = kk+1;\r
		end\r
		t = 0.1^(kk);\r
		pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi));\r
\r
		%Pour phi et rho\r
		for r=1:k\r
			for mm=1:m\r
				for i=1:n\r
					ps1(i,mm,r) = Y2(i,mm,r) * dot(X2(i,:,r), phi(:,mm,r));\r
					nY21(i,mm,r) = (Y2(i,mm,r))^2;\r
				end\r
				ps(mm,r) = sum(ps1(:,mm,r));\r
				nY2(mm,r) = sum(nY21(:,mm,r));\r
				rho(mm,mm,r) = ((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));\r
			end\r
		end\r
		for r=1:k\r
			for j=1:p\r
				for mm=1:m	\r
					S(j,mm,r) = -rho(mm,mm,r)*ps2(j,mm,r) + dot(phi(1:j-1,mm,r),Gram2(j,1:j-1,r)')...\r
						+ dot(phi(j+1:p,mm,r),Gram2(j,j+1:p,r)');\r
					if abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma)\r
						phi(j,mm,r)=0;\r
					else \r
						if S(j,mm,r)> n*lambda*(pi(r)^gamma)\r
							phi(j,mm,r)=(n*lambda*(pi(r)^gamma)-S(j,mm,r))/Gram2(j,j,r);\r
						else\r
							phi(j,mm,r)=-(n*lambda*(pi(r)^gamma)+S(j,mm,r))/Gram2(j,j,r);\r
						end\r
					end\r
				end\r
			end\r
		end\r
\r
		%%%%%%%%%%\r
		%Etape E %\r
		%%%%%%%%%%\r
		\r
		sumLogLLF2 = 0.0;\r
		for i=1:n\r
			%precompute dot products to numerically adjust their values\r
			dotProducts = zeros(k,1);\r
			for r=1:k\r
				dotProducts(r)= (Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r));\r
			end\r
			shift = 0.5*min(dotProducts);\r
			\r
			%compute Gam(:,:) using shift determined above\r
			sumLLF1 = 0.0;\r
			for r=1:k\r
				Gam(i,r) = pi(r)*det(rho(:,:,r))*exp(-0.5*dotProducts(r) + shift);\r
				sumLLF1 = sumLLF1 + Gam(i,r)/(2*PI)^(m/2);\r
			end\r
			sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
			sumGamI = sum(Gam(i,:));\r
			if sumGamI > EPS\r
				gam(i,:) = Gam(i,:) / sumGamI;\r
			else\r
				gam(i,:) = zeros(k,1);\r
			end\r
		end\r
		\r
		sumPen = 0.0;\r
		for r=1:k\r
			sumPen = sumPen + pi(r).^gamma .* b(r);\r
		end\r
		LLF(ite) = -(1/n)*sumLogLLF2 + lambda*sumPen;\r
		\r
		if ite == 1\r
			dist = LLF(ite);\r
		else \r
			dist = (LLF(ite)-LLF(ite-1))/(1+abs(LLF(ite)));\r
		end\r
		\r
		Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))));\r
		Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))));\r
		Dist3=max(max((abs(pi-Pi))./(1+abs(Pi))));\r
		dist2=max([Dist1,Dist2,Dist3]); \r
		\r
		ite=ite+1;\r
	end\r
\r
	pi = transpose(pi);\r
\r
end\r