src/test/generate_test_data/helpers/constructionModelesLassoMLE.R

   1 constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2){
   2   #get matrix sizes
   3   n = dim(X)[1];
   4   p = dim(phiInit)[1]
   5   m = dim(phiInit)[2]
   6   k  = dim(phiInit)[3]
   7   L = length(glambda)
   8
   9   #output parameters
  10   phi = array(0, dim=c(p,m,k,L))
  11   rho = array(0, dim=c(m,m,k,L))
  12   Pi = matrix(0, k, L)
  13   llh = matrix(0, L, 2) #log-likelihood
  14
  15   for(lambdaIndex in 1:L){
  16     a = A1[, 1, lambdaIndex]
  17     a[a==0] = c()
  18     if(length(a)==0){
  19       next
  20     }
  21     EMGLLf = EMGLLF(phiInit[a,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0,X[,a],Y,tau)
  22
  23     phiLambda = EMGLLf$phi
  24     rhoLambda = EMGLLf$rho
  25     piLambda = EMGLLf$Pi
  26
  27     for(j in 1:length(a)){
  28       phi[a[j],,,lambdaIndex] = phiLambda[j,,]
  29     }
  30     rho[,,,lambdaIndex] = rhoLambda
  31     Pi[,lambdaIndex] = piLambda
  32
  33     dimension = 0
  34     for(j in 1:p){
  35       vec =  c(2, dim(A2)[2])
  36       b = A2[j,vec,lambdaIndex]
  37       b[b==0] = c()
  38       if(length(b) > 0){
  39         phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0
  40       }
  41       c = A1[j,vec,lambdaIndex]
  42       c[c==0] = c()
  43       dimension = dimension + length(c)
  44     }
  45
  46     #on veut calculer l'EMV avec toutes nos estimations
  47         densite = matrix(0, n, L)
  48         for(i in 1:n){
  49             for( r in 1:k){
  50                 delta = Y[i,]%*%rho[,,r,lambdaIndex] - (X[i,a]%*%phi[a,,r,lambdaIndex]);
  51                 densite[i,lambdaIndex] = densite[i,lambdaIndex] +   Pi[r,lambdaIndex]*det(rho[,,r,lambdaIndex])/(sqrt(2*pi))^m*exp(-tcrossprod(delta)/2.0)
  52             }
  53         }
  54         llh[lambdaIndex,1] = sum(log(densite[,lambdaIndex]))
  55         llh[lambdaIndex,2] = (dimension+m+1)*k-1
  56   }
  57   return(list(phi=phi, rho=rho, Pi=Pi, llh = llh))
  58 }