-library(MASS) #generalized inverse of matrix Monroe-Penrose
-
-vec_bin = function(X,r){
- Z = c()
- indice = c()
- j=1
- for(i in 1:length(X)){
- if(X[i] == r){
- Z[i] = 1
- indice[j] = i
- j=j+1
- }
- else{
- Z[i] = 0
- }
- }
- return(list(Z,indice))
+vec_bin = function(X,r)
+{
+ Z = c()
+ indice = c()
+ j = 1
+ for (i in 1:length(X))
+ {
+ if(X[i] == r)
+ {
+ Z[i] = 1
+ indice[j] = i
+ j=j+1
+ } else
+ Z[i] = 0
+ }
+ return (list(Z=Z,indice=indice))
}
-initSmallEM = function(k,X,Y,tau){
- n = nrow(Y)
- m = ncol(Y)
- p = ncol(X)
+initSmallEM = function(k,X,Y,tau)
+{
+ n = nrow(Y)
+ m = ncol(Y)
+ p = ncol(X)
- betaInit1 = array(0, dim=c(p,m,k,20))
- sigmaInit1 = array(0, dim = c(m,m,k,20))
- phiInit1 = array(0, dim = c(p,m,k,20))
- rhoInit1 = array(0, dim = c(m,m,k,20))
- piInit1 = matrix(0,20,k)
- gamInit1 = array(0, dim=c(n,k,20))
- LLFinit1 = list()
-
-
- for(repet in 1:20){
- clusters = hclust(dist(y)) #default distance : euclidean
- clusterCut = cutree(clusters,k)
- Zinit1[,repet] = clusterCut #retourne les indices (à quel cluster indiv_i appartient) d'un clustering hierarchique (nb de cluster = k)
-
- for(r in 1:k){
- Z = Zinit1[,repet]
- Z_bin = vec_bin(Z,r)
- Z_vec = Z_bin[[1]] #vecteur 0 et 1 aux endroits où Z==r
- Z_indice = Z_bin[[2]] #renvoit les indices où Z==r
-
- betaInit1[,,r,repet] = ginv(t(x[Z_indice,])%*%x[Z_indice,])%*%t(x[Z_indice,])%*%y[Z_indice,]
- sigmaInit1[,,r,repet] = diag(m)
- phiInit1[,,r,repet] = betaInit1[,,r,repet]/sigmaInit1[,,r,repet]
- rhoInit1[,,r,repet] = solve(sigmaInit1[,,r,repet])
- piInit1[repet,r] = sum(Z_vec)/n
- }
-
- for(i in 1:n){
- for(r in 1:k){
- dotProduct = (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet]) %*% (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet])
- Gam[i,r] = piInit1[repet,r]*det(rhoInit1[,,r,repet])*exp(-0.5*dotProduct)
- }
- sumGamI = sum(gam[i,])
- gamInit1[i,,repet]= Gam[i,] / sumGamI
- }
-
- miniInit = 10
- maxiInit = 11
-
- new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,x,y,tau)
- ##.C("EMGLLF", phiInit = phiInit, rhoInit = rhoInit, ...)
- LLFEessai = new_EMG[[4]]
- LLFinit1[[repet]] = LLFEessai[[length(LLFEessai)]]
- }
-
- b = which.max(LLFinit1)
-
- phiInit = phiInit1[,,,b]
- rhoInit = rhoInit1[,,,b]
- piInit = piInit1[b,]
- gamInit = gamInit1[,,b]
-
- return(list(phiInit, rhoInit, piInit, gamInit))
-}
+ betaInit1 = array(0, dim=c(p,m,k,20))
+ sigmaInit1 = array(0, dim = c(m,m,k,20))
+ phiInit1 = array(0, dim = c(p,m,k,20))
+ rhoInit1 = array(0, dim = c(m,m,k,20))
+ piInit1 = matrix(0,20,k)
+ gamInit1 = array(0, dim=c(n,k,20))
+ LLFinit1 = list()
+
+ require(MASS) #Moore-Penrose generalized inverse of matrix
+ for(repet in 1:20)
+ {
+ clusters = hclust(dist(y)) #default distance : euclidean
+ #cutree retourne les indices (à quel cluster indiv_i appartient) d'un clustering hierarchique
+ clusterCut = cutree(clusters,k)
+ Zinit1[,repet] = clusterCut
+
+ for(r in 1:k)
+ {
+ Z = Zinit1[,repet]
+ Z_bin = vec_bin(Z,r)
+ Z_vec = Z_bin$Z #vecteur 0 et 1 aux endroits où Z==r
+ Z_indice = Z_bin$indice #renvoit les indices où Z==r
+
+ betaInit1[,,r,repet] =
+ ginv(t(x[Z_indice,])%*%x[Z_indice,])%*%t(x[Z_indice,])%*%y[Z_indice,]
+ sigmaInit1[,,r,repet] = diag(m)
+ phiInit1[,,r,repet] = betaInit1[,,r,repet]/sigmaInit1[,,r,repet]
+ rhoInit1[,,r,repet] = solve(sigmaInit1[,,r,repet])
+ piInit1[repet,r] = sum(Z_vec)/n
+ }
+ for(i in 1:n)
+ {
+ for(r in 1:k)
+ {
+ dotProduct = (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet]) %*%
+ (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet])
+ Gam[i,r] = piInit1[repet,r]*det(rhoInit1[,,r,repet])*exp(-0.5*dotProduct)
+ }
+ sumGamI = sum(gam[i,])
+ gamInit1[i,,repet]= Gam[i,] / sumGamI
+ }
+ miniInit = 10
+ maxiInit = 11
+
+ new_EMG = .Call("EMGLLF",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],
+ gamInit1[,,repet],miniInit,maxiInit,1,0,x,y,tau)
+ LLFEessai = new_EMG$LLF
+ LLFinit1[repet] = LLFEessai[length(LLFEessai)]
+ }
+
+ b = which.max(LLFinit1)
+ phiInit = phiInit1[,,,b]
+ rhoInit = rhoInit1[,,,b]
+ piInit = piInit1[b,]
+ gamInit = gamInit1[,,b]
+
+ return (list(phiInit=phiInit, rhoInit=rhoInit, piInit=piInit, gamInit=gamInit))
+}
-## 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')
+SelMix = setRefClass(
+ Class = "selmix",
+
+ fields = c(
+ # User defined
+
+ # regression data (size n*p, where n is the number of observations,
+ # and p is the number of regressors)
+ X = "numeric",
+ # response data (size n*m, where n is the number of observations,
+ # and m is the number of responses)
+ Y = "numeric",
+
+ # Optionally user defined (some default values)
+
+ # power in the penalty
+ gamma = "double",
+ # minimum number of iterations for EM algorithm
+ mini = "integer",
+ # maximum number of iterations for EM algorithm
+ maxi = "integer",
+ # threshold for stopping EM algorithm
+ eps = "double",
+ # minimum number of components in the mixture
+ kmin = "integer",
+ # maximum number of components in the mixture
+ kmax = "integer",
+ rangmin = "integer",
+ rangmax = "integer",
+
+ # Computed through the workflow
+
+ # initialisation for the reparametrized conditional mean parameter
+ phiInit,
+ # initialisation for the reparametrized variance parameter
+ rhoInit,
+ # initialisation for the proportions
+ piInit,
+ # initialisation for the allocations probabilities in each component
+ tauInit,
+ # values for the regularization parameter grid
+ gridLambda = [];
+ # 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 ...
+ A1,
+ A2,
+ # collection of estimations for the reparametrized conditional mean parameters
+ Phi,
+ # collection of estimations for the reparametrized variance parameters
+ Rho,
+ # collection of estimations for the proportions parameters
+ Pi,
+
+ #immutable
+ seuil = 1e-15;
+ ),
+
+ methods = list(
+ #######################
+ #initialize main object
+ #######################
+ initialize = function(X,Y,...)
+ {
+ "Initialize SelMix object"
+
+ callSuper(...)
+
+ X <<- X;
+ Y <<- Y;
+ gamma <<- ifelse (hasArg("gamma"), gamma, 1.)
+ mini <<- ifelse (hasArg("mini"), mini, as.integer(5))
+ maxi <<- ifelse (hasArg("maxi"), maxi, as.integer(10))
+ eps <<- ifelse (hasArg("eps"), eps, 1e-6)
+ kmin <<- ifelse (hasArg("kmin"), kmin, as.integer(2))
+ kmax <<- ifelse (hasArg("kmax"), kmax, as.integer(3))
+ rangmin <<- ifelse (hasArg("rangmin"), rangmin, as.integer(2))
+ rangmax <<- ifelse (hasArg("rangmax"), rangmax, as.integer(3))
+ },
+
+ ##################################
+ #core workflow: compute all models
+ ##################################
+
+ initParameters = function(k)
+ {
+ "Parameters initialization"
+
+ #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.
+ init = initSmallEM(k,sx.X,sx.Y,sx.eps);
+ phiInit <<- init$phi0;
+ rhoInit <<- init$rho0;
+ piInit <<- init$pi0;
+ tauInit <<- init$tau0;
+ },
+
+ computeGridLambda = function()
+ {
+ "computation of the regularization grid"
+ #(according to explicit formula given by EM algorithm)
+
+ gridLambda <<- grillelambda(sx.phiInit,sx.rhoInit,sx.piInit,sx.tauInit,sx.X,sx.Y,
+ sx.gamma,sx.mini,sx.maxi,sx.eps);
+ },
+
+ computeRelevantParameters = function()
+ {
+ "Compute relevant parameters"
+
+ #select variables according to each regularization parameter
+ #from the grid: sx.A1 corresponding to selected variables, and
+ #sx.A2 corresponding to unselected variables.
+ params = 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);
+ A1 <<- params$A1
+ A2 <<- params$A2
+ Rho <<- params$Rho
+ Pi <<- params$Pi
+ },
+
+ runProcedure1 = function()
+ {
+ "Run procedure 1 [EMGLLF]"
+
+ #compute parameter estimations, with the Maximum Likelihood
+ #Estimator, restricted on selected variables.
+ res = 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);
+ return (list( phi=res$phi, rho=res$rho, pi=res$pi))
+ },
+
+ runProcedure2 = function()
+ {
+ "Run procedure 2 [EMGrank]"
+
+ #compute parameter estimations, with the Low Rank
+ #Estimator, restricted on selected variables.
+ return (constructionModelesLassoRank(sx.Pi,sx.Rho,sx.mini,sx.maxi,sx.X,sx.Y,sx.eps,
+ sx.A1,sx.rangmin,sx.rangmax)$phi)
+ },
+
+ run = function(procedure)
+ {
+ "main loop: over all k and all lambda"
+
+ # Run the all procedure, 1 with the
+ #maximum likelihood refitting, and 2 with the Low Rank refitting.
+ p = dim(phiInit)[1]
+ m = dim(phiInit)[2]
+ for (k in kmin:kmax)
+ {
+ print(k)
+ initParameters(k)
+ computeGridLambda()
+ computeRelevantParameters()
+ if (procedure == 1)
+ {
+ r1 = runProcedure1(sx)
+ Phi2 = Phi
+ Rho2 = Rho
+ Pi2 = Pi
+ p = ncol(X)
+ m = ncol(Y)
+ if size(Phi2) == 0 #TODO: continue translation MATLAB --> R
+ sx.Phi(:,:,1:k,:) = r1$phi;
+ sx.Rho(:,:,1:k,:) = r1$rho;
+ sx.Pi(1:k,:) = r1$pi;
+ else
+ sx.Phi = zeros(p,m,sx.kmax,size(Phi2,4)+size(r1$phi,4));
+ sx.Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2;
+ sx.Phi(:,:,1:k,size(Phi2,4)+1:end) = r1$phi;
+ sx.Rho = zeros(m,m,sx.kmax,size(Rho2,4)+size(r1$rho,4));
+ sx.Rho(:,:,1:size(Rho2,3),1:size(Rho2,4)) = Rho2;
+ sx.Rho(:,:,1:k,size(Rho2,4)+1:end) = r1$rho;
+ sx.Pi = zeros(sx.kmax,size(Pi2,2)+size(r1$pi,2));
+ sx.Pi(1:size(Pi2,1),1:size(Pi2,2)) = Pi2;
+ sx.Pi(1:k,size(Pi2,2)+1:end) = r1$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
+
+ )
+)
projects / valse.git / commitdiff