X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=src%2Ftest%2Fgenerate_test_data%2Fhelpers%2FEMGLLF.R;h=7c80081e09fc18a9beadff926428f6a406262a9c;hp=02fbb47e7476f25cb245af27656219dc4ef34881;hb=ef67d338c7f28ba041abe40ca9a8ab69f8365a90;hpb=c3bc47052f3ccb659659c59a82e9a99ea842398d diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R index 02fbb47..7c80081 100644 --- a/src/test/generate_test_data/helpers/EMGLLF.R +++ b/src/test/generate_test_data/helpers/EMGLLF.R @@ -1,4 +1,5 @@ -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] @@ -8,11 +9,10 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) #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)) @@ -23,7 +23,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) 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] ##bon calcul ? idem pour X2 ??... - } - 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){ + for (j in 1:p) + { + for (s in 1:p) Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r]) - } } } @@ -58,116 +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]))) - } + for (r in 1:k) + b[r] = sum(abs(phi[,,r])) gam2 = colSums(gam) - a = sum(gam%*%(log(Pi))) + a = sum(gam %*% log(pi)) #tant que les props sont negatives kk = 0 pi2AllPositive = FALSE - while(pi2AllPositive == FALSE){ - Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi) - pi2AllPositive = TRUE - for(r in 1:k){ - if(Pi2[r] < 0){ - pi2AllPositive = false; - break - } - } + while (!pi2AllPositive) + { + pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi) + 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)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(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] * (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])) + 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){ - p1 = p-1 - for(j in 1:p1){ - for(mm in 1:m){ - j1 = j-1 - j2 = j+1 - v1 = c(1:j1) - v2 = c(j2:p) - S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r]) + phi[v2,mm,r]%*%(Gram2[j,v2,r]) #erreur indice - 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{ - 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] } } } - + ########## #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(c(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 = t(Pi) - return(list("phi"=phi, "rho"=rho, "pi"=Pi, "LLF"=LLF, "S"=S)) + + return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S)) }