From c202886908a6eff2cd704f910a0bde06be5a0875 Mon Sep 17 00:00:00 2001 From: Benjamin Goehry Date: Tue, 24 Jan 2017 18:07:07 +0100 Subject: [PATCH] EMGrank + typo --- src/test/generate_test_data/helpers/EMGLLF.R | 27 +++---- src/test/generate_test_data/helpers/EMGrank.R | 72 +++++++++++++++++++ .../generate_test_data/helpers/checkOutput.R | 6 +- .../generate_test_data/helpers/covariance.R | 4 +- 4 files changed, 90 insertions(+), 19 deletions(-) create mode 100644 src/test/generate_test_data/helpers/EMGrank.R diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R index f108a38..94a917d 100644 --- a/src/test/generate_test_data/helpers/EMGLLF.R +++ b/src/test/generate_test_data/helpers/EMGLLF.R @@ -1,9 +1,9 @@ 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 @@ -68,10 +68,10 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) kk = 0 pi2AllPositive = FALSE while(pi2AllPositive == FALSE){ - pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi) + Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi) pi2AllPositive = TRUE for(r in 1:k){ - if(pi2[r] < 0){ + if(Pi2[r] < 0){ pi2AllPositive = false; break } @@ -81,12 +81,12 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) #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) + 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 = 0.1^(kk) - pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi)) + Pi = (Pi+t*(Pi2-Pi)) / sum(Pi+t*(Pi2-Pi)) #Pour phi et rho for(r in 1:k){ @@ -104,13 +104,14 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) 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)) + if(abs(S(j,mm,r)) <= n*lambda*(Pi[r]^gamma)){ phi[j,mm,r]=0 - else{ - if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)) + }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] + }else{ + phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r] + } } } } diff --git a/src/test/generate_test_data/helpers/EMGrank.R b/src/test/generate_test_data/helpers/EMGrank.R new file mode 100644 index 0000000..0ee2d1f --- /dev/null +++ b/src/test/generate_test_data/helpers/EMGrank.R @@ -0,0 +1,72 @@ +EMGLLF = function(Pi, Rho, mini, maxi, X, Y, tau, rank){ + #matrix dimensions + n = dim(X)[1] + p = dim(X)[2] + m = dim(Rho)[2] + k = dim(Rho)[3] + + #init outputs + phi = array(0, dim=c(p,m,k)) + Z = rep(1, n) + Pi = piInit + LLF = 0 + + #local variables + Phi = array(0, dim=c(p,m,k)) + deltaPhi = c(0) + sumDeltaPhi = 0 + deltaPhiBufferSize = 20 + + #main loop + ite = 1 + while(ite<=mini || (ite<=maxi && sumDeltaPhi>tau)){ + #M step: Mise à jour de Beta (et donc phi) + for(r in 1:k){ + Z_bin = vec_bin(Z,r) + Z_vec = Z_bin$vec #vecteur 0 et 1 aux endroits o? Z==r + Z_indice = Z_bin$indice + if(sum(Z_indice) == 0){ + next + } + #U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr + [U,S,V] = svd(ginv(crossprod(X[Z_indice,]))%*% (X[Z_indice,])%*%Y[Z_indice,] ) + #Set m-rank(r) singular values to zero, and recompose + #best rank(r) approximation of the initial product + S[rank(r)+1:end,] = 0 + phi[,,r] = U %*%S%*%t(V)%*%Rho[,,r] + } + + #Etape E et calcul de LLF + sumLogLLF2 = 0 + for(i in 1:n){ + sumLLF1 = 0 + maxLogGamIR = -Inf + for(r in 1:k){ + dotProduct = tcrossprod(Y[i,]%*%Rho[,,r]-X[i,]%*%phi[,,r]) + logGamIR = log(Pi[r]) + log(det(Rho[,,r])) - 0.5*dotProduct + #Z[i] = index of max (gam[i,]) + if(logGamIR > maxLogGamIR){ + Z[i] = r + maxLogGamIR = logGamIR + } + sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2) + } + sumLogLLF2 = sumLogLLF2 + log(sumLLF1) + } + + LLF = -1/n * sumLogLLF2 + + #update distance parameter to check algorithm convergence (delta(phi, Phi)) + deltaPhi = c(deltaPhi, max(max(max((abs(phi-Phi))/(1+abs(phi))))) ) + if(length(deltaPhi) > deltaPhiBufferSize){ + deltaPhi = deltaPhi[2:length(deltaPhi)] + } + sumDeltaPhi = sum(abs(deltaPhi)) + + #update other local variables + Phi = phi + ite = ite+1 + + } + return(list(phi=phi, LLF=LLF)) +} \ No newline at end of file diff --git a/src/test/generate_test_data/helpers/checkOutput.R b/src/test/generate_test_data/helpers/checkOutput.R index 187536a..ae2fbdf 100644 --- a/src/test/generate_test_data/helpers/checkOutput.R +++ b/src/test/generate_test_data/helpers/checkOutput.R @@ -2,11 +2,9 @@ checkOutput = function(varName, array, refArray, tol) { print(paste("Checking ",varName,sep="")) maxError = max(abs(array - refArray)) - if(maxError >= tol) - { + if(maxError >= tol){ print(paste("Inaccuracy: max(abs(error)) = ",maxError," >= ",tol,sep="")) - } else - { + } else{ print("OK") } } diff --git a/src/test/generate_test_data/helpers/covariance.R b/src/test/generate_test_data/helpers/covariance.R index 15cd693..09a9ec5 100644 --- a/src/test/generate_test_data/helpers/covariance.R +++ b/src/test/generate_test_data/helpers/covariance.R @@ -1,8 +1,8 @@ covariance = function(p,a) { A = matrix(a, p,p) - for (i in 1:p) + for (i in 1:p){ A[i,] = A[i,]^abs(i-(1:p)) - + } return (A) } -- 2.44.0