R/constructionModelesLassoMLE.R

   1 constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,
   2     X,Y,seuil,tau,A1,A2)
   3 {
   4   n = dim(X)[1];
   5   p = dim(phiInit)[1]
   6   m = dim(phiInit)[2]
   7   k = dim(phiInit)[3]
   8   L = length(glambda)
   9
  10   #output parameters
  11   phi = array(0, dim=c(p,m,k,L))
  12   rho = array(0, dim=c(m,m,k,L))
  13   pi = matrix(0, k, L)
  14   llh = matrix(0, L, 2) #log-likelihood
  15
  16   for(lambdaIndex in 1:L)
  17     {
  18     a = A1[,1,lambdaIndex]
  19     a = a[a!=0]
  20     if(length(a)==0)
  21       next
  22
  23     res = EMGLLF(phiInit[a,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0.,X[,a],Y,tau)
  24
  25         #TODO: supprimer ça et utiliser parLapply(...)
  26     for (j in 1:length(a))
  27       phi[a[j],,,lambdaIndex] = res$phi[j,,]
  28     rho[,,,lambdaIndex] = res$rho
  29     pi[,lambdaIndex] = res$pi
  30
  31     dimension = 0
  32     for (j in 1:p) #TODO: doit pouvoir être fait beaucoup plus simplement
  33         {
  34       b = A2[j,2:dim(A2)[2],lambdaIndex]
  35       b = b[b!=0]
  36       if (length(b) > 0)
  37         phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0.
  38       c = A1[j,2:dim(A1)[2],lambdaIndex]
  39       dimension = dimension + sum(c!=0)
  40     }
  41
  42     #on veut calculer l'EMV avec toutes nos estimations
  43         densite = matrix(0, nrow=n, ncol=L)
  44         for (i in 1:n) #TODO: pas besoin de cette boucle (vectorize)
  45         {
  46             for (r in 1:k)
  47             {
  48                 delta = Y[i,]%*%rho[,,r,lambdaIndex] - (X[i,a]%*%phi[a,,r,lambdaIndex]);
  49                 densite[i,lambdaIndex] = densite[i,lambdaIndex] + pi[r,lambdaIndex] *
  50                     det(rho[,,r,lambdaIndex])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0)
  51             }
  52         }
  53         llh[lambdaIndex,1] = sum(log(densite[,lambdaIndex]))
  54         llh[lambdaIndex,2] = (dimension+m+1)*k-1
  55   }
  56   return (list("phi"=phi, "rho"=rho, "pi"=pi, "llh" = llh))
  57 }