X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=src%2Ftest%2Fgenerate_test_data%2Fhelpers%2FEMGLLF.R;h=7c80081e09fc18a9beadff926428f6a406262a9c;hb=ef67d338c7f28ba041abe40ca9a8ab69f8365a90;hp=f108a38e7d27d478cc55b7af8ef781e33d4ff10c;hpb=83ed2c0a1924fe874d0d12f82ad46c42b6ae4e97;p=valse.git diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R index f108a38..7c80081 100644 --- a/src/test/generate_test_data/helpers/EMGLLF.R +++ b/src/test/generate_test_data/helpers/EMGLLF.R @@ -1,29 +1,29 @@ -EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){ +EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) +{ #matrix dimensions n = dim(X)[1] - p = dim[phiInit][1] - m = dim[phiInit][2] - k = dim[phiInit][3] + p = dim(phiInit)[1] + m = dim(phiInit)[2] + k = dim(phiInit)[3] #init outputs phi = phiInit rho = rhoInit - Pi = piInit + pi = piInit LLF = rep(0, maxi) S = array(0, dim=c(p,m,k)) - gam = gamInit Gram2 = array(0, dim=c(p,p,k)) ps2 = array(0, dim=c(p,m,k)) b = rep(0, k) pen = matrix(0, maxi, k) X2 = array(0, dim=c(n,p,k)) - Y2 = array(0, dim=c(p,m,k)) + Y2 = array(0, dim=c(n,m,k)) dist = 0 dist2 = 0 ite = 1 - Pi2 = rep(0, k) + pi2 = rep(0, k) ps = matrix(0, m,k) nY2 = matrix(0, m,k) ps1 = array(0, dim=c(n,m,k)) @@ -31,25 +31,25 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) Gam = matrix(0, n,k) EPS = 1E-15 - while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))){ + while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))) + { Phi = phi Rho = rho - PI = Pi + Pi = pi + #calcul associé à Y et X - for(r in 1:k){ - for(mm in 1:m){ - Y2[,mm,r] = sqrt(gam[,r]) .* Y[,mm] - } - for(i in 1:n){ - X2[i,,r] = X[i,] .* sqrt(gam[i,r]) - } - for(mm in 1:m){ + for(r in 1:k) + { + for (mm in 1:m) + Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm] + for (i in 1:n) + X2[i,,r] = sqrt(gam[i,r]) * X[i,] + for (mm in 1:m) ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r]) - } - for(j in 1:p){ - for(s in 1:p){ - Gram2[j,s,r] = tcrossprod(X2[,j,r], X2[,s,r]) - } + for (j in 1:p) + { + for (s in 1:p) + Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r]) } } @@ -58,110 +58,106 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) ########## #pour pi - for(r in 1:k){ - b[r] = sum(sum(abs(phi[,,r]))) - } - gam2 = sum(gam[1,]) #BIG DOUTE - a = sum(gam*t(log(Pi))) + for (r in 1:k) + b[r] = sum(abs(phi[,,r])) + gam2 = colSums(gam) + a = sum(gam %*% log(pi)) #tant que les props sont negatives kk = 0 pi2AllPositive = FALSE - while(pi2AllPositive == FALSE){ + while (!pi2AllPositive) + { pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi) - pi2AllPositive = TRUE - for(r in 1:k){ - if(pi2[r] < 0){ - pi2AllPositive = false; - break - } - } + pi2AllPositive = all(pi2 >= 0) kk = kk+1 } - #t[m]la plus grande valeur dans la grille O.1^k tel que ce soit - #décroissante ou constante - while((-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*t(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000){ - pi2 = pi+0.1^kk*(1/n*gam2-pi) - kk = kk+1 + #t[m] la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante + while( kk < 1000 && -a/n + lambda * sum(pi^gamma * b) < + -sum(gam2 * log(pi2))/n + lambda * sum(pi2^gamma * b) ) + { + pi2 = pi + 0.1^kk * (1/n*gam2 - pi) + kk = kk + 1 } - t = 0.1^(kk) - pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi)) + t = 0.1^kk + pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi)) #Pour phi et rho - for(r in 1:k){ - for(mm in 1:m){ - for(i in 1:n){ - ps1[i,mm,r] = Y2[i,mm,r] * dot(X2(i,:,r), phi(:,mm,r)) - nY21[i,mm,r] = (Y2[i,mm,r])^2 + for (r in 1:k) + { + for (mm in 1:m) + { + for (i in 1:n) + { + ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r]) + nY21[i,mm,r] = Y2[i,mm,r]^2 } - ps[mm,r] = sum(ps1(:,mm,r)); - nY2[mm,r] = sum(nY21(:,mm,r)); - rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r])) + ps[mm,r] = sum(ps1[,mm,r]) + nY2[mm,r] = sum(nY21[,mm,r]) + rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r]) } } - for(r in 1:k){ - for(j in 1:p){ - for(mm in 1:m){ - 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]) + dot(phi[j+1:p,mm,r],Gram2[j,j+1:p,r]) - if(abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma)) + for (r in 1:k) + { + for (j in 1:p) + { + for (mm in 1:m) + { + S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + + (if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) + + (if(j n*lambda*(Pi[r]^gamma)) - phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r] + else if(S[j,mm,r] > n*lambda*(pi[r]^gamma)) + phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r] else - phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r] - } + phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r] } } } - + ########## #Etape E # ########## sumLogLLF2 = 0 - for(i in 1:n){ - #precompute dot products to numerically adjust their values - dotProducts = rep(0,k) - for(r in 1:k){ - dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r]) - } - shift = 0.5*min(dotProducts) - + for (i in 1:n) + { + #precompute sq norms to numerically adjust their values + sqNorm2 = rep(0,k) + for (r in 1:k) + sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 ) + shift = 0.5*min(sqNorm2) + #compute Gam(:,:) using shift determined above sumLLF1 = 0.0; - for(r in 1:k){ - Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift) - sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2) + for (r in 1:k) + { + #FIXME: numerical problems, because 0 < det(Rho[,,r] < EPS; what to do ?! + # consequence: error in while() at line 77 + Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r] + shift) * det(rho[,,r]) + sumLLF1 = sumLLF1 + Gam[i,r] / (2*base::pi)^(m/2) } sumLogLLF2 = sumLogLLF2 + log(sumLLF1) sumGamI = sum(Gam[i,]) if(sumGamI > EPS) gam[i,] = Gam[i,] / sumGamI else - gam[i,] = rep(0,k) - } - - - sumPen = 0 - for(r in 1:k){ - sumPen = sumPen + Pi[r].^gamma^b[r] + gam[i,] = rep(0,k) } - LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen - - if(ite == 1) - dist = LLF[ite] - else - dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite])) - - Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi))))) - Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho))))) - Dist3=max(max((abs(Pi-PI))./(1+abs(PI)))) - dist2=max([Dist1,Dist2,Dist3]) - - ite=ite+1 + + sumPen = sum(pi^gamma * b) + LLF[ite] = -sumLogLLF2/n + lambda*sumPen + + dist = ifelse( ite == 1, LLF[ite], (LLF[ite]-LLF[ite-1]) / (1+abs(LLF[ite])) ) + + Dist1 = max( (abs(phi-Phi)) / (1+abs(phi)) ) + Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) ) + Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)) ) + dist2 = max(Dist1,Dist2,Dist3) + + ite = ite+1 } - - Pi = transpose(Pi) - return(list(phi=phi, rho=rho, Pi=Pi, LLF=LLF, S=S)) -} \ No newline at end of file + + return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S)) +}