## TODO: turn this code into R classdef selmix < handle properties (SetAccess = private) %always user defined X; % regression data (size n*p, where n is the number of observations, and p is the number of regressors) Y; % response data (size n*m, where n is the number of observations, and m is the number of responses) %optionally user defined some default values gamma; % power in the penalty mini; % minimum number of iterations for EM algorithm maxi; % maximum number of iterations for EM algorithm eps; % threshold for stopping EM algorithm kmin; % minimum number of components in the mixture kmax; % maximum number of components in the mixture rangmin; rangmax; %computed through the workflow phiInit; % initialisation for the reparametrized conditional mean parameter rhoInit; % initialisation for the reparametrized variance parameter piInit; % initialisation for the proportions tauInit; % initialisation for the allocations probabilities in each component gridLambda = []; % values for the regularization parameter grid 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 ... A2; Phi; % collection of estimations for the reparametrized conditional mean parameters Rho; % collection of estimations for the reparametrized variance parameters Pi; % collection of estimations for the proportions parameters end properties (Constant) %immutable seuil = 1e-15; end methods %%%%%%%%%%%%%%%%%%%%%%% %initialize main object %%%%%%%%%%%%%%%%%%%%%%% function sx = selmix(X,Y,varargin) %set defaults for optional inputs optargs = {1.0 5 10 1e-6 2 3 2 3}; %replace defaults by user parameters optargs(1:length(varargin)) = varargin; sx.X = X; sx.Y = Y; [sx.gamma,sx.mini,sx.maxi,sx.eps,sx.kmin,sx.kmax,sx.rangmin,sx.rangmax] = optargs{:}; % z = somme(sx.X,sx.Y); % sx.Z=z; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %core workflow: compute all models %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function initParameters(sx,k) [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. sx.phiInit = phi0; sx.rhoInit = rho0; sx.piInit = pi0; sx.tauInit = tau0; end function computeGridLambda(sx) [sx.gridLambda] = grillelambda(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.X,sx.Y,sx.gamma,sx.mini,sx.maxi,sx.eps); % computation of the regularization grid, according to explicit % formulae given by EM algorithm. end function computeRelevantParameters(sx) [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); %select variables according to each regularization parameter %from the grid: sx.A1 corresponding to selected variables, and %sx.A2 corresponding to unselected variables. end function [sx,phi,rho,pi]=runProcedure1(sx) [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); %compute parameter estimations, with the Maximum Likelihood %Estimator, restricted on selected variables. end function [phi] =runProcedure2(sx) [phi,~] = constructionModelesLassoRank(sx.Pi,sx.Rho,sx.mini,sx.maxi,sx.X,sx.Y,sx.eps,sx.A1,sx.rangmin,sx.rangmax); %compute parameter estimations, with the Low Rank %Estimator, restricted on selected variables. end % main loop: over all k and all lambda function run(sx,procedure) % Run the all procedure, 1 with the %maximum likelihood refitting, and 2 with the Low Rank refitting. [p,m,~]=size(sx.phiInit); for k=sx.kmin:sx.kmax k initParameters(sx,k); computeGridLambda(sx); computeRelevantParameters(sx); if (procedure == 1) [~,phi,rho,pi] = runProcedure1(sx); Phi2 = sx.Phi; Rho2 = sx.Rho; Pi2 = sx.Pi; p = size(sx.X,2); m = size(sx.Y,2); if size(Phi2) == 0 sx.Phi(:,:,1:k,:) = phi; sx.Rho(:,:,1:k,:) = rho; sx.Pi(1:k,:) = pi; else sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4)); sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2; sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi; sx.Rho = zeros(m,m,sx.kmax,size(Rho2,4)+size(rho,4)); sx.Rho(:,:,1:size(Rho2,3),1:size(Rho2,4)) = Rho2; sx.Rho(:,:,1:k,size(Rho2,4)+1:end) = rho; sx.Pi = zeros(sx.kmax,size(Pi2,2)+size(pi,2)); sx.Pi(1:size(Pi2,1),1:size(Pi2,2)) = Pi2; sx.Pi(1:k,size(Pi2,2)+1:end) = pi; end else [phi] = runProcedure2(sx); phi Phi2 = sx.Phi; if size(Phi2,1) == 0 sx.Phi(:,:,1:k,:) = phi; else size(Phi2) sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(phi,4)); size(sx.Phi) sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2; sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = phi; end end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %pruning: select only one (or a few best ?!) model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % function[model] selectModel(sx) % %TODO % %model = sxModel(...); % end end end %%%%%%%%%%%%% %OLD VERSION: %~ for k=K %~ % On initialise %~ initialisation %~ disp('Initialisé') %~ % On construit la grille des lambdas : variables informatives %~ [glambda,phiEmv,rhoEmv,piEmv]=grillelambda(phiInit,rhoInit,piInit,tauInit,x,y,gamma,mini,maxi,tau); %~ glambda=glambda(1:3); %~ disp('glambda construite') %~ % On trouve les variables informatives pour chaque lambda : S est la %~ % matrice des coefficients des variables informatives %~ % on parallelise à l interieur du selectiontotale() %~ [B1,B2,rhoLasso,piLasso]=selectiontotale(phiInit,rhoInit,piInit,tauInit,mini,maxi,gamma,glambda,x,y,10^-15,tau); %~ %S1 les variables informatives, S2 celles qui ne le sont pas %~ [B1bis,B2bis,glambda2bis,ind,rhoLasso,piLasso]=suppressionmodelesegaux(B1,B2,glambda,rhoLasso,piLasso); %~ dessinVariablesSelectionnees; %~ disp('Le Lasso est fait') %~ % Pour chaque lambda ainsi construit, on veut calculer l'EMV pour la procédure Lasso-MLE %~ %On obtient une collection de modèles pour Lasso-MLE %~ % ICI AUSSI ON PEUT PARALLELISER a l interieur de constructionModelesLassoMLE %~ nombreLambda=size(B1bis,3); %~ %[phiLassoMLE,rhoLassoMLE,piLassoMLE,vraisLassoMLE]=constructionModelesLassoMLE(phiInit,rhoInit,piInit, tauInit,mini,maxi,gamma,glambda2bis,X,Y,seuil,tau,B1bis,B2bis) %~ %Pour Lasso-Rank %~ %on ne garde que les colonnes actives %~ B1ter=B1bis(:,1,:); %~ % [B1ter,ind,rhoLasso,piLasso]=suppressionmodelesegaux2(B1ter,rhoLasso,piLasso) %~ %Pour chaque lambda, on veut calculer l'estimateur low rank pour la procédure Lasso-Rank %~ %On obtient une collection de modèles pour Lasso-Rank %~ %ICI AUSSI ON PEUT PARALLELISER a linterieur de constructionModelesLassoRank %~ nombreLambda2=size(B1ter,2); %~ [phi,lvraisemblance,Z] = constructionModelesLassoRank(... %~ piEmv,rhoEmv,mini,maxi,X,Y,tau,B1ter,2,4); %~ disp('On a construit les modèles pour les deux procédures') %~ end %~ %selectionModelesLassoMLE; %~ %selectionModelesLassoRank; %~ disp('On a sélectionné les modèles dans les deux procédures')