R/main.R

   1 ## TODO: turn this code into R
   2
   3 classdef selmix < handle
   4
   5     properties (SetAccess = private)
   6         %always user defined
   7         X; % regression data (size n*p, where n is the number of observations, and p is the number of regressors)
   8         Y; % response data (size n*m, where n is the number of observations, and m is the number of responses)
   9
  10         %optionally user defined some default values
  11         gamma; % power in the penalty
  12         mini; % minimum number of iterations for EM algorithm
  13         maxi; % maximum number of iterations for EM algorithm
  14         eps; % threshold for stopping EM algorithm
  15         kmin; % minimum number of components in the mixture
  16         kmax; % maximum number of components in the mixture
  17         rangmin;
  18         rangmax;
  19
  20         %computed through the workflow
  21         phiInit; % initialisation for the reparametrized conditional mean parameter
  22         rhoInit; % initialisation for the reparametrized variance parameter
  23         piInit; % initialisation for the proportions
  24         tauInit; % initialisation for the allocations probabilities in each component
  25         gridLambda = []; % values for the regularization parameter grid
  26         A1; %je ne crois pas vraiment qu'il faille les mettre en sortie, d'autant plus qu'on construit une matrice A1 et A2 pour chaque k, et elles sont grandes, donc ca coute un peu cher ...
  27         A2;
  28         Phi; % collection of estimations for the reparametrized conditional mean parameters
  29         Rho; % collection of estimations for the reparametrized variance parameters
  30         Pi; % collection of estimations for the proportions parameters
  31     end
  32
  33     properties (Constant)
  34         %immutable
  35         seuil = 1e-15;
  36     end
  37
  38     methods
  39         %%%%%%%%%%%%%%%%%%%%%%%
  40         %initialize main object
  41         %%%%%%%%%%%%%%%%%%%%%%%
  42
  43         function sx = selmix(X,Y,varargin)
  44             %set defaults for optional inputs
  45             optargs = {1.0 5 10 1e-6 2 3 2 3};
  46             %replace defaults by user parameters
  47             optargs(1:length(varargin)) = varargin;
  48             sx.X = X;
  49             sx.Y = Y;
  50             [sx.gamma,sx.mini,sx.maxi,sx.eps,sx.kmin,sx.kmax,sx.rangmin,sx.rangmax] = optargs{:};
  51             %           z = somme(sx.X,sx.Y);
  52             %           sx.Z=z;
  53         end
  54
  55
  56         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  57         %core workflow: compute all models
  58         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  59
  60         function initParameters(sx,k)
  61             [phi0,rho0,pi0,tau0] = initSmallEM(k,sx.X,sx.Y,sx.eps); %smallEM initializes parameters by k-means and regression model in each component, doing this 20 times, and keeping the values maximizing the likelihood after 10 iterations of the EM algorithm.
  62             sx.phiInit = phi0;
  63             sx.rhoInit = rho0;
  64             sx.piInit  = pi0;
  65             sx.tauInit = tau0;
  66         end
  67
  68         function computeGridLambda(sx)
  69             [sx.gridLambda] = grillelambda(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.X,sx.Y,sx.gamma,sx.mini,sx.maxi,sx.eps);
  70             % computation of the regularization grid, according to explicit
  71             % formulae given by EM algorithm.
  72         end
  73
  74         function computeRelevantParameters(sx)
  75             [sx.A1,sx.A2,sx.Rho,sx.Pi] = selectiontotale(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.mini,sx.maxi,sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.eps);
  76             %select variables according to each regularization parameter
  77             %from the grid: sx.A1 corresponding to selected variables, and
  78             %sx.A2 corresponding to unselected variables.
  79         end
  80
  81         function [sx,phi,rho,pi]=runProcedure1(sx)
  82             [phi,rho,pi,~] = constructionModelesLassoMLE(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.mini,sx.maxi,sx.gamma,sx.gridLambda,sx.X,sx.Y,sx.seuil,sx.eps,sx.A1,sx.A2);
  83             %compute parameter estimations, with the Maximum Likelihood
  84             %Estimator, restricted on selected variables.
  85         end
  86
  87         function [phi] =runProcedure2(sx)
  88             [phi,~] = constructionModelesLassoRank(sx.Pi,sx.Rho,sx.mini,sx.maxi,sx.X,sx.Y,sx.eps,sx.A1,sx.rangmin,sx.rangmax);
  89             %compute parameter estimations, with the Low Rank
  90             %Estimator, restricted on selected variables.
  91         end
  92
  93         % main loop: over all k and all lambda
  94         function run(sx,procedure) % Run the all procedure, 1 with the
  95             %maximum likelihood refitting, and 2 with the Low Rank refitting.
  96             [p,m,~]=size(sx.phiInit);
  97             for k=sx.kmin:sx.kmax
  98                 k
  99                 initParameters(sx,k);
 100                 computeGridLambda(sx);
 101                 computeRelevantParameters(sx);
 102                 if (procedure == 1)
 103                     [~,phi,rho,pi] = runProcedure1(sx);
 104                     Phi2 = sx.Phi;
 105                     Rho2 = sx.Rho;
 106                     Pi2 = sx.Pi;
 107                     p = size(sx.X,2);
 108                     m = size(sx.Y,2);
 109                     if size(Phi2) == 0
 110                         sx.Phi(:,:,1:k,:) = phi;
 111                         sx.Rho(:,:,1:k,:) = rho;
 112                         sx.Pi(1:k,:) = pi;
 113                     else
 114                         sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4));
 115                         sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
 116                         sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi;
 117                         sx.Rho = zeros(m,m,sx.kmax,size(Rho2,4)+size(rho,4));
 118                         sx.Rho(:,:,1:size(Rho2,3),1:size(Rho2,4)) = Rho2;
 119                         sx.Rho(:,:,1:k,size(Rho2,4)+1:end) = rho;
 120                         sx.Pi = zeros(sx.kmax,size(Pi2,2)+size(pi,2));
 121                         sx.Pi(1:size(Pi2,1),1:size(Pi2,2)) = Pi2;
 122                         sx.Pi(1:k,size(Pi2,2)+1:end) = pi;
 123                     end
 124                 else
 125                     [phi] = runProcedure2(sx);
 126                     phi
 127                     Phi2 = sx.Phi;
 128                     if size(Phi2,1) == 0
 129                         sx.Phi(:,:,1:k,:) = phi;
 130                     else
 131                         size(Phi2)
 132                         sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4));
 133                         size(sx.Phi)
 134                         sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
 135                         sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi;
 136                     end
 137
 138                 end
 139
 140
 141             end
 142         end
 143
 144         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 145         %pruning: select only one (or a few best ?!) model
 146         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 147         %
 148         %       function[model] selectModel(sx)
 149         %           %TODO
 150         %           %model = sxModel(...);
 151         %       end
 152
 153     end
 154
 155 end
 156
 157
 158 %%%%%%%%%%%%%
 159 %OLD VERSION:
 160
 161 %~ for k=K
 162 %~ % On initialise
 163 %~ initialisation
 164 %~ disp('Initialisé')
 165 %~ % On construit la grille des lambdas : variables informatives
 166 %~ [glambda,phiEmv,rhoEmv,piEmv]=grillelambda(phiInit,rhoInit,piInit,tauInit,x,y,gamma,mini,maxi,tau);
 167 %~ glambda=glambda(1:3);
 168 %~ disp('glambda construite')
 169 %~ % On trouve les variables informatives pour chaque lambda : S est la
 170 %~ % matrice des coefficients des variables informatives
 171 %~ % on parallelise à l interieur du selectiontotale()
 172 %~ [B1,B2,rhoLasso,piLasso]=selectiontotale(phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,glambda,x,y,10^-15,tau);
 173 %~ %S1 les variables informatives, S2 celles qui ne le sont pas
 174 %~ [B1bis,B2bis,glambda2bis,ind,rhoLasso,piLasso]=suppressionmodelesegaux(B1,B2,glambda,rhoLasso,piLasso);
 175 %~ dessinVariablesSelectionnees;
 176 %~ disp('Le Lasso est fait')
 177 %~ % Pour chaque lambda ainsi construit, on veut calculer l'EMV pour la procédure Lasso-MLE
 178 %~ %On obtient une collection de modèles pour Lasso-MLE
 179 %~ % ICI AUSSI ON PEUT PARALLELISER a l interieur de constructionModelesLassoMLE
 180 %~ nombreLambda=size(B1bis,3);
 181 %~ %[phiLassoMLE,rhoLassoMLE,piLassoMLE,vraisLassoMLE]=constructionModelesLassoMLE(phiInit,rhoInit,piInit, tauInit,mini,maxi,gamma,glambda2bis,X,Y,seuil,tau,B1bis,B2bis)
 182 %~ %Pour Lasso-Rank
 183 %~ %on ne garde que les colonnes actives
 184 %~ B1ter=B1bis(:,1,:);
 185 %~ %        [B1ter,ind,rhoLasso,piLasso]=suppressionmodelesegaux2(B1ter,rhoLasso,piLasso)
 186 %~ %Pour chaque lambda, on veut calculer l'estimateur low rank pour la procédure Lasso-Rank
 187 %~ %On obtient une collection de modèles pour Lasso-Rank
 188 %~ %ICI AUSSI ON PEUT PARALLELISER a linterieur de constructionModelesLassoRank
 189 %~ nombreLambda2=size(B1ter,2);
 190 %~ [phi,lvraisemblance,Z] = constructionModelesLassoRank(...
 191 %~ piEmv,rhoEmv,mini,maxi,X,Y,tau,B1ter,2,4);
 192 %~ disp('On a construit les modèles pour les deux procédures')
 193 %~ end
 194 %~ %selectionModelesLassoMLE;
 195 %~ %selectionModelesLassoRank;
 196 %~ disp('On a sélectionné les modèles dans les deux procédures')