[valse.git] / R / main.R

## 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')
Commit	Line	Data
493a35bf BA	1	## TODO: turn this code into R
493a35bf BA	2
1d3c1faa BA	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')