+++ /dev/null
-function[phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n,p,m,k)
-
- phiInit = zeros(p,m,k);
-
- piInit = (1.0/k) * ones(1,k);
-
- rhoInit = zeros(m,m,k);
- for r=1:k
- rhoInit(:,:,r) = eye(m,m);
- end
-
- gamInit = 0.1 * ones(n,k);
- R = random('unid',k,n,1);
- for i=1:n
- gamInit(i,R(i)) = 0.9;
- end
- gamInit = gamInit / (sum(gamInit(1,:)));
-
-end
+++ /dev/null
-%compile C code (for MATLAB or Octave)
-if exist('octave_config_info')
- setenv('CFLAGS','-O2 -std=gnu99 -fopenmp')
- mkoctfile --mex -DOctave -I../Util -I../ProcLassoMLE selectiontotale.c selectiontotale_interface.c ../Util/ioutils.c ../ProcLassoMLE/EMGLLF.c -o selectiontotale -lm -lgsl -lgslcblas -lgomp
-else
- mex CFLAGS="\$CFLAGS -O2 -std=gnu99 -fopenmp" -I../Util -I../ProcLassoMLE selectiontotale.c selectiontotale_interface.c ../Util/ioutils.c ../ProcLassoMLE/EMGLLF.c -output selectiontotale -lm -lgsl -lgslcblas -lgomp
-end
+++ /dev/null
-%X is generated following a gaussian mixture \sum pi_r N(meanX_k, covX_r)
-%Y is generated then, with Y_i ~ \sum pi_r N(Beta_r.X_i, covY_r)
-function[X,Y,Z] = generateIO(meanX, covX, covY, pi, beta, n)
-
- [p, ~, k] = size(covX);
- [m, ~, ~] = size(covY);
- if exist('octave_config_info')
- %Octave statistics package doesn't have gmdistribution()
- X = zeros(n, p);
- Z = zeros(n);
- cs = cumsum(pi);
- for i=1:n
- %TODO: vectorize ? http://stackoverflow.com/questions/2977497/weighted-random-numbers-in-matlab
- tmpRand01 = rand();
- [~,Z(i)] = min(cs - tmpRand01 >= 0);
- X(i,:) = mvnrnd(meanX(Z(i),:), covX(:,:,Z(i)), 1);
- end
- else
- gmDistX = gmdistribution(meanX, covX, pi);
- [X, Z] = random(gmDistX, n);
- end
-
- Y = zeros(n, m);
- BX = zeros(n,m,k);
- for i=1:n
- for r=1:k
- %compute beta_r . X_i
- BXir = zeros(1, m);
- for mm=1:m
- BXir(mm) = dot(X(i,:), beta(:,mm,r));
- end
- %add pi(r) * N(beta_r . X_i, covY) to Y_i
- Y(i,:) = Y(i,:) + pi(r) * mvnrnd(BXir, covY(:,:,r), 1);
- end
- end
-
-end
+++ /dev/null
-%call generateIO with default parameters (random means, covariances = identity, equirepartition)
-function[X,Y,Z] = generateIOdefault(n, p, m, k)
-
- rangeX = 100;
- meanX = rangeX * (1 - 2*rand(k, p));
- covX = zeros(p,p,k);
- covY = zeros(m,m,k);
- for r=1:k
- covX(:,:,r) = eye(p);
- covY(:,:,r) = eye(m);
- end
- pi = (1/k) * ones(1,k);
-
- %initialize beta to a random number of non-zero random value
- beta = zeros(p,m,k);
- for j=1:p
- nonZeroCount = ceil(m*rand(1));
- beta(j,1:nonZeroCount,:) = rand(nonZeroCount, k);
- end
-
- [X,Y,Z] = generateIO(meanX, covX, covY, pi, beta, n);
-
-end
+++ /dev/null
-function[gridLambda] = grillelambda(phiInit,rhoInit,piInit,gamInit,X,Y,gamma,mini,maxi,tau)
-
- n = size(X, 1);
- [p,m,k] = size(phiInit);
- [phi,rho,pi,~,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,1,0,X,Y,tau);
-
- gridLambda = zeros(p,m,k);
- for j=1:p
- for mm=1:m
- gridLambda(j,mm,:) = abs(reshape(S(j,mm,:),[k,1])) ./ (n*pi.^gamma);
- end
- end
-
- gridLambda = unique(gridLambda);
- gridLambda(gridLambda()>1) = [];
-
-end
+++ /dev/null
-function[phiInit,rhoInit,piInit,gamInit] = initSmallEM(k,x,y,tau)
-[n,m]=size(y);
-gamInit1=zeros(n,k,20);
-for repet=1:20
- Zinit1(:,repet)=clusterdata(y,k);
- for r=1:k
- betaInit1(:,:,r,repet)=pinv(transpose(x(Zinit1(:,repet)==r,:))*x(Zinit1(:,repet)==r,:))*transpose(x(Zinit1(:,repet)==r,:))*y(Zinit1(:,repet)==r,:);
- sigmaInit1(:,:,r,repet)=eye(m);
- phiInit1(:,:,r,repet)=betaInit1(:,:,r,repet)/sigmaInit1(:,:,r,repet);
- rhoInit1(:,:,r,repet)=inv(sigmaInit1(:,:,r,repet));
- piInit1(repet,r)=sum(Zinit1(:,repet)==r)/n;
- end
- for i=1:n
- for r=1:k
- dotProduct = (y(i,:)*rhoInit1(:,:,r,repet)-x(i,:)*phiInit1(:,:,r,repet)) * transpose(y(i,:)*rhoInit1(:,:,r,repet)-x(i,:)*phiInit1(:,:,r,repet));
- Gam(i,r) = piInit1(repet,r)*det(rhoInit1(:,:,r,repet))*exp(-0.5*dotProduct);
- end
- sumGamI = sum(Gam(i,:));
- gamInit1(i,:,repet) = Gam(i,:) / sumGamI;
- end
- miniInit=int64(10);
- maxiInit=int64(11);
-
- [~,~,~,LLFEssai,~] = EMGLLF(phiInit1(:,:,:,repet),rhoInit1(:,:,:,repet),piInit1(repet,:),gamInit1(:,:,repet),miniInit,maxiInit,1,0,x,y,tau);
- LLFinit1(repet)=LLFEssai(end);
-end
-[~,b]=max(LLFinit1);
-
-phiInit=phiInit1(:,:,:,b);
-rhoInit=rhoInit1(:,:,:,b);
-piInit=piInit1(b,:);
-gamInit=gamInit1(:,:,b);
-end
\ No newline at end of file
+++ /dev/null
-function[A1,A2,Rho,Pi] = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau)
-
- [p,m,k] = size(phiInit);
- L = length(glambda);
- A1 = zeros(p,m+1,L,'int64');
- A2 = zeros(p,m+1,L,'int64');
- Rho = zeros(m,m,k,L);
- Pi = zeros(k,L);
-
- %Pour chaque lambda de la grille, on calcule les coefficients
- for lambdaIndex=1:L
- [phi,rho,pi,~,~] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda(lambdaIndex),X,Y,tau);
-
- %Si un des coefficients est supérieur au seuil, on garde cette variable
- selectedVariables = zeros(p,m);
- discardedVariables = zeros(p,m);
- atLeastOneSelectedVariable = false;
- for j=1:p
- cpt=1;
- cpt2=1;
- for mm=1:m
- if max(abs(phi(j,mm,:))) > seuil
- selectedVariables(j,cpt) = mm;
- cpt = cpt+1;
- atLeastOneSelectedVariable = true;
- else
- discardedVariables(j,cpt2) = mm;
- cpt2 = cpt2+1;
- end
- end
- end
-
- %Si aucun des coefficients n'a été gardé on renvoit la matrice nulle
- %Et si on enlevait ces colonnes de zéro ??? Indices des colonnes vides
- if atLeastOneSelectedVariable
- vec = [];
- for j=1:p
- if selectedVariables(j,1) ~= 0
- vec = [vec;j];
- end
- end
-
- %Sinon on renvoit les numéros des coefficients utiles
- A1(:,1,lambdaIndex) = [vec;zeros(p-length(vec),1)];
- A1(1:length(vec),2:m+1,lambdaIndex) = selectedVariables(vec,:);
- A2(:,1,lambdaIndex) = 1:p;
- A2(:,2:m+1,lambdaIndex) = discardedVariables;
- Rho(:,:,:,lambdaIndex) = rho;
- Pi(:,lambdaIndex) = pi;
- end
-
- end
-
-end
+++ /dev/null
-function[phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)\r
-\r
- %Get matrices dimensions\r
- PI = 4.0 * atan(1.0);\r
- n = size(X, 1);\r
- [p,m,k] = size(phiInit);\r
- \r
- %Initialize outputs\r
- phi = phiInit;\r
- rho = rhoInit;\r
- pi = piInit;\r
- LLF = zeros(maxi,1);\r
- S = zeros(p,m,k);\r
- \r
- %Other local variables\r
- %NOTE: variables order is always n,p,m,k\r
- gam = gamInit;\r
- Gram2 = zeros(p,p,k);\r
- ps2 = zeros(p,m,k);\r
- b = zeros(k,1);\r
- pen = zeros(maxi,k);\r
- X2 = zeros(n,p,k);\r
- Y2 = zeros(n,m,k);\r
- dist = 0;\r
- dist2 = 0;\r
- ite = 1;\r
- pi2 = zeros(k,1);\r
- ps = zeros(m,k);\r
- nY2 = zeros(m,k);\r
- ps1 = zeros(n,m,k);\r
- nY21 = zeros(n,m,k);\r
- Gam = zeros(n,k);\r
- EPS = 1e-15;\r
-\r
- while ite<=mini || (ite<=maxi && (dist>=tau || dist2>=sqrt(tau)))\r
- \r
- Phi = phi;\r
- Rho = rho;\r
- Pi = pi;\r
- \r
- %Calculs associés à Y et X\r
- for r=1:k\r
- for mm=1:m\r
- Y2(:,mm,r) = sqrt(gam(:,r)) .* Y(:,mm);\r
- end\r
- for i=1:n\r
- X2(i,:,r) = X(i,:) .* sqrt(gam(i,r));\r
- end\r
- for mm=1:m\r
- ps2(:,mm,r) = transpose(X2(:,:,r)) * Y2(:,mm,r);\r
- end\r
- for j=1:p\r
- for s=1:p\r
- Gram2(j,s,r) = dot(X2(:,j,r), X2(:,s,r));\r
- end\r
- end\r
- end\r
-\r
- %%%%%%%%%%\r
- %Etape M %\r
- %%%%%%%%%%\r
- \r
- %Pour pi\r
- for r=1:k\r
- b(r) = sum(sum(abs(phi(:,:,r))));\r
- end\r
- gam2 = sum(gam,1);\r
- a = sum(gam*transpose(log(pi)));\r
- \r
- %tant que les proportions sont negatives\r
- kk = 0;\r
- pi2AllPositive = false;\r
- while ~pi2AllPositive\r
- pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi);\r
- pi2AllPositive = true;\r
- for r=1:k\r
- if pi2(r) < 0\r
- pi2AllPositive = false;\r
- break;\r
- end\r
- end\r
- kk = kk+1;\r
- end\r
- \r
- %t(m) la plus grande valeur dans la grille O.1^k tel que ce soit\r
- %décroissante ou constante\r
- while (-1/n*a+lambda*((pi.^gamma)*b))<(-1/n*gam2*transpose(log(pi2))+lambda.*(pi2.^gamma)*b) && kk<1000\r
- pi2 = pi+0.1^kk*(1/n*gam2-pi);\r
- kk = kk+1;\r
- end\r
- t = 0.1^(kk);\r
- pi = (pi+t*(pi2-pi)) / sum(pi+t*(pi2-pi));\r
-\r
- %Pour phi et rho\r
- for r=1:k\r
- for mm=1:m\r
- for i=1:n\r
- ps1(i,mm,r) = Y2(i,mm,r) * dot(X2(i,:,r), phi(:,mm,r));\r
- nY21(i,mm,r) = (Y2(i,mm,r))^2;\r
- end\r
- ps(mm,r) = sum(ps1(:,mm,r));\r
- nY2(mm,r) = sum(nY21(:,mm,r));\r
- rho(mm,mm,r) = ((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));\r
- end\r
- end\r
- for r=1:k\r
- for j=1:p\r
- for mm=1:m \r
- S(j,mm,r) = -rho(mm,mm,r)*ps2(j,mm,r) + dot(phi(1:j-1,mm,r),Gram2(j,1:j-1,r)')...\r
- + dot(phi(j+1:p,mm,r),Gram2(j,j+1:p,r)');\r
- if abs(S(j,mm,r)) <= n*lambda*(pi(r)^gamma)\r
- phi(j,mm,r)=0;\r
- else \r
- if S(j,mm,r)> n*lambda*(pi(r)^gamma)\r
- phi(j,mm,r)=(n*lambda*(pi(r)^gamma)-S(j,mm,r))/Gram2(j,j,r);\r
- else\r
- phi(j,mm,r)=-(n*lambda*(pi(r)^gamma)+S(j,mm,r))/Gram2(j,j,r);\r
- end\r
- end\r
- end\r
- end\r
- end\r
-\r
- %%%%%%%%%%\r
- %Etape E %\r
- %%%%%%%%%%\r
- \r
- sumLogLLF2 = 0.0;\r
- for i=1:n\r
- %precompute dot products to numerically adjust their values\r
- dotProducts = zeros(k,1);\r
- for r=1:k\r
- dotProducts(r)= (Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*rho(:,:,r)-X(i,:)*phi(:,:,r));\r
- end\r
- shift = 0.5*min(dotProducts);\r
- \r
- %compute Gam(:,:) using shift determined above\r
- sumLLF1 = 0.0;\r
- for r=1:k\r
- Gam(i,r) = pi(r)*det(rho(:,:,r))*exp(-0.5*dotProducts(r) + shift);\r
- sumLLF1 = sumLLF1 + Gam(i,r)/(2*PI)^(m/2);\r
- end\r
- sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
- sumGamI = sum(Gam(i,:));\r
- if sumGamI > EPS\r
- gam(i,:) = Gam(i,:) / sumGamI;\r
- else\r
- gam(i,:) = zeros(k,1);\r
- end\r
- end\r
- \r
- sumPen = 0.0;\r
- for r=1:k\r
- sumPen = sumPen + pi(r).^gamma .* b(r);\r
- end\r
- LLF(ite) = -(1/n)*sumLogLLF2 + lambda*sumPen;\r
- \r
- if ite == 1\r
- dist = LLF(ite);\r
- else \r
- dist = (LLF(ite)-LLF(ite-1))/(1+abs(LLF(ite)));\r
- end\r
- \r
- Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))));\r
- Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))));\r
- Dist3=max(max((abs(pi-Pi))./(1+abs(Pi))));\r
- dist2=max([Dist1,Dist2,Dist3]); \r
- \r
- ite=ite+1;\r
- end\r
-\r
- pi = transpose(pi);\r
-\r
-end\r
+++ /dev/null
-%compile C code (for MATLAB or Octave)
-if exist('octave_config_info')
- setenv('CFLAGS','-O2 -std=gnu99 -fopenmp')
- mkoctfile --mex -DOctave -I../Util EMGLLF.c EMGLLF_interface.c ../Util/ioutils.c -o EMGLLF -lm -lgsl -lgslcblas -lgomp
- mkoctfile --mex -DOctave -I../Util constructionModelesLassoMLE.c constructionModelesLassoMLE_interface.c EMGLLF.c ../Util/ioutils.c -o constructionModelesLassoMLE -lm -lgsl -lgslcblas -lgomp
-else
- mex CFLAGS="\$CFLAGS -O2 -std=gnu99 -fopenmp" -I../Util EMGLLF.c EMGLLF_interface.c ../Util/ioutils.c -output EMGLLF -lm -lgsl -lgslcblas -lgomp
- mex CFLAGS="\$CFLAGS -O2 -std=gnu99 -fopenmp" -I../Util constructionModelesLassoMLE.c constructionModelesLassoMLE_interface.c EMGLLF.c ../Util/ioutils.c -output constructionModelesLassoMLE -lm -lgsl -lgslcblas -lgomp
-end
+++ /dev/null
-function[phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(...
- phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2)
-
- PI = 4.0 * atan(1.0);
-
- %get matrix sizes
- n = size(X, 1);
- [p,m,k] = size(phiInit);
- L = length(glambda);
-
- %output parameters
- phi = zeros(p,m,k,L);
- rho = zeros(m,m,k,L);
- pi = zeros(k,L);
- lvraisemblance = zeros(L,2);
-
- for lambdaIndex=1:L
- % Procedure Lasso-MLE
- a = A1(:,1,lambdaIndex);
- a(a==0) = [];
- if length(a) == 0
- continue;
- end
- [phiLambda,rhoLambda,piLambda,~,~] = EMGLLF(...
- phiInit(a,:,:),rhoInit,piInit,gamInit,mini,maxi,gamma,0,X(:,a),Y,tau);
-
- for j=1:length(a)
- phi(a(j),:,:,lambdaIndex) = phiLambda(j,:,:);
- end
- rho(:,:,:,lambdaIndex) = rhoLambda;
- pi(:,lambdaIndex) = piLambda;
-
- dimension = 0;
- for j=1:p
- b = A2(j,2:end,lambdaIndex);
- b(b==0) = [];
- if length(b) > 0
- phi(A2(j,1,lambdaIndex),b,:,lambdaIndex) = 0.0;
- end
- c = A1(j,2:end,lambdaIndex);
- c(c==0) = [];
- dimension = dimension + length(c);
- end
-
- %on veut calculer l'EMV avec toutes nos estimations
- densite = zeros(n,L);
- for i=1:n
- for r=1:k
- delta = Y(i,:)*rho(:,:,r,lambdaIndex) - (X(i,a)*(phi(a,:,r,lambdaIndex)));
- densite(i,lambdaIndex) = densite(i,lambdaIndex) +...
- pi(r,lambdaIndex)*det(rho(:,:,r,lambdaIndex))/(sqrt(2*PI))^m*exp(-dot(delta,delta)/2.0);
- end
- end
- lvraisemblance(lambdaIndex,1) = sum(log(densite(:,lambdaIndex)));
- lvraisemblance(lambdaIndex,2) = (dimension+m+1)*k-1;
- end
-
-end
+++ /dev/null
-function[phi,LLF] = EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank)\r
-\r
- % get matrix sizes\r
- [~,m,k] = size(Rho);\r
- [n,p] = size(X);\r
-\r
- % allocate output matrices\r
- phi = zeros(p,m,k);\r
- Z = ones(n,1,'int64');\r
- LLF = 0.0;\r
-\r
- % local variables\r
- Phi = zeros(p,m,k);\r
- deltaPhi = 0.0;\r
- deltaPhi = [];\r
- sumDeltaPhi = 0.0;\r
- deltaPhiBufferSize = 20;\r
-\r
- %main loop (at least mini iterations)\r
- ite = int64(1);\r
- while ite<=mini || (ite<=maxi && sumDeltaPhi>tau)\r
-\r
- %M step: Mise à jour de Beta (et donc phi)\r
- for r=1:k\r
- if (sum(Z==r) == 0)\r
- continue;\r
- end\r
- %U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr\r
- [U,S,V] = svd(pinv(transpose(X(Z==r,:))*X(Z==r,:))*transpose(X(Z==r,:))*Y(Z==r,:));\r
- %Set m-rank(r) singular values to zero, and recompose \r
- %best rank(r) approximation of the initial product\r
- S(rank(r)+1:end,:) = 0;\r
- phi(:,:,r) = U * S * transpose(V) * Rho(:,:,r);\r
- end\r
- \r
- %Etape E et calcul de LLF\r
- sumLogLLF2 = 0.0;\r
- for i=1:n\r
- sumLLF1 = 0.0;\r
- maxLogGamIR = -Inf;\r
- for r=1:k\r
- dotProduct = (Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r)) * transpose(Y(i,:)*Rho(:,:,r)-X(i,:)*phi(:,:,r));\r
- logGamIR = log(Pi(r)) + log(det(Rho(:,:,r))) - 0.5*dotProduct;\r
- %Z(i) = index of max (gam(i,:))\r
- if logGamIR > maxLogGamIR\r
- Z(i) = r;\r
- maxLogGamIR = logGamIR;\r
- end\r
- sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2);\r
- end\r
- sumLogLLF2 = sumLogLLF2 + log(sumLLF1);\r
- end\r
- \r
- LLF = -1/n * sumLogLLF2;\r
-\r
- % update distance parameter to check algorithm convergence (delta(phi, Phi))\r
- deltaPhi = [ deltaPhi, max(max(max((abs(phi-Phi))./(1+abs(phi))))) ];\r
- if length(deltaPhi) > deltaPhiBufferSize\r
- deltaPhi = deltaPhi(2:length(deltaPhi));\r
- end\r
- sumDeltaPhi = sum(abs(deltaPhi));\r
-\r
- % update other local variables\r
- Phi = phi;\r
- ite = ite+1;\r
-\r
- end\r
-\r
-end\r
+++ /dev/null
-%compile C code (for MATLAB or Octave)
-if exist('octave_config_info')
- setenv('CFLAGS','-O2 -std=gnu99 -fopenmp')
- mkoctfile --mex -DOctave -I../Util EMGrank.c EMGrank_interface.c ../Util/ioutils.c -o EMGrank -lm -lgsl -lgslcblas -lgomp
- mkoctfile --mex -DOctave -I../Util constructionModelesLassoRank.c constructionModelesLassoRank_interface.c EMGrank.c ../Util/ioutils.c -o constructionModelesLassoRank -lm -lgsl -lgslcblas -lgomp
-else
- mex CFLAGS="\$CFLAGS -O2 -std=gnu99 -fopenmp" -I../Util EMGrank.c EMGrank_interface.c ../Util/ioutils.c -output EMGrank -lm -lgsl -lgslcblas -lgomp
- mex CFLAGS="\$CFLAGS -O2 -std=gnu99 -fopenmp" -I../Util constructionModelesLassoRank.c constructionModelesLassoRank_interface.c EMGrank.c ../Util/ioutils.c -output constructionModelesLassoRank -lm -lgsl -lgslcblas -lgomp
-end
+++ /dev/null
-function[phi,lvraisemblance] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
-
- PI = 4.0 * atan(1.0);
-
- %get matrix sizes
- [n,p] = size(X);
- [~,m,k,~] = size(Rho);
- L = size(A1, 2); %A1 est p x m+1 x L ou p x L ?!
-
- %On cherche les rangs possiblement intéressants
- deltaRank = rangmax - rangmin + 1;
- Size = deltaRank^k;
- Rank = zeros(Size,k,'int64');
- for r=1:k
- %On veut le tableau de toutes les combinaisons de rangs possibles
- %Dans la première colonne : on répète (rangmax-rangmin)^(k-1) chaque chiffre : ca remplit la colonne
- %Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre, et on fait ca (rangmax-rangmin)^2 fois
- %...
- %Dans la dernière, on répète chaque chiffre une fois, et on fait ca (rangmin-rangmax)^(k-1) fois.
- Rank(:,r) = rangmin + reshape(repmat(0:(deltaRank-1), deltaRank^(k-r), deltaRank^(r-1)), Size, 1);
- end
-
- %output parameters
- phi = zeros(p,m,k,L*Size);
- lvraisemblance = zeros(L*Size,2);
- for lambdaIndex=1:L
- %On ne garde que les colonnes actives
- %active sera l'ensemble des variables informatives
- active = A1(:,lambdaIndex);
- active(active==0) = [];
- if length(active) > 0
- for j=1:Size
- [phiLambda,LLF] = EMGrank(Pi(:,lambdaIndex),Rho(:,:,:,lambdaIndex),mini,maxi,X(:,active),Y,tau,Rank(j,:));
- lvraisemblance((lambdaIndex-1)*Size+j,:) = [LLF, sum(Rank(j,:) .* (length(active)-Rank(j,:)+m))];
- phi(active,:,:,(lambdaIndex-1)*Size+j) = phiLambda;
- end
- end
- end
-
-end
+++ /dev/null
-%Selection de modèle dans la procédure Lasso-MLE
-% Selection de modele : voici les parametres selectionnes pour chaque
-% dimension, et l'indice associe
-[indiceLassoMLE,D1LassoMLE] = selectionmodele(vraisLassoMLE);
-
-%on veut tracer la courbe
-%figure(2)
-%plot(D1LassoMLE/n,vraisLassoMLE(indiceLassoMLE),'.')
-%on veut appliquer Capushe : il nous faut un tableau de données
-tableauLassoMLE = enregistrerdonnees(D1LassoMLE,vraisLassoMLE(indiceLassoMLE,:),n);
-save tableauLassoMLE.mat tableauLassoMLE
-%On conclut avec Capushe !
+++ /dev/null
-%Selection de modèle dans la procedure Lasso Rank
-vraisLassoRank=vraisLassoRank';
-% Selection de modele : voici les parametres selectionnes pour chaque
-% dimension, et l'indice associe
-[indiceLassoRank,D1LassoRank]=selectionmodele(vraisLassoRank);
-tableauLassoRank=enregistrerdonnees(D1LassoRank,vraisLassoRank(indiceLassoRank,:),n);
-%On veut tracer la courbe des pentes
-%figure(2)
-%plot(D1LassoRank/n,vraisLassoRank(indiceLassoRank),'.')
-save tableauLassoRank.mat tableauLassoRank
-%On conclut avec Capushe !
+++ /dev/null
-%utile dans selectiontotale.m, remarque quels coefficients sont nuls et lesquels ne le sont pas
-function[A,B]=selectiondindice(phi,seuil)
-
- [pp,m,~]=size(phi);
- A=zeros(pp,m);
- B=zeros(pp,m);
- for j=1:pp
- cpt=0;cpt2=0;
- for mm=1:m
- if (max(phi(j,mm,:)) > seuil)
- cpt=cpt+1;
- A(j,cpt)=mm;
- else cpt2=cpt2+1;
- B(j,cpt2)=mm;
- end
- end
- end
-
-end
+++ /dev/null
-function [indice,D1]=selectionmodele(vraisemblance)
-
- D=vraisemblance(:,2);
- [D1]=unique(D);
- indice=ones(1,length(D1));
- %On ne sélectionne que celui qui maximise : l'EMV
- if length(D1)>2
- for i=1:length(D1)
- a=[];
- for j=1:length(D)
- if D(j)==D1(i)
- a=[a,vraisemblance(j,1)];
- end
- end
- b=max(a);
- indice(i)=find(vraisemblance(:,1)==b,1);
- end
- end
-
-end
+++ /dev/null
-function[B1,B2,glambda,ind,rho,pi]=suppressionmodelesegaux(B1,B2,glambda,rho,pi)
-
- ind=[];
- for l=1:length(glambda)
- for ll=1:l-1
- if B1(:,:,l)==B1(:,:,ll)
- ind=[ind l];
- end
- end
- end
- ind=unique(ind);
- B1(:,:,ind)=[];
- glambda(ind)=[];
- B2(:,:,ind)=[];
- rho(:,:,:,ind)=[];
- pi(:,ind)=[];
-
-end
+++ /dev/null
-function[B1,ind,rho,pi]=suppressionmodelesegaux2(B1,rho,pi)
-
- ind=[];
- nombreLambda=size(B1,2);
- for l=1:nombreLambda
- for ll=1:l-1
- if B1(:,l)==B1(:,ll)
- ind=[ind l];
- end
- end
- end
- ind=unique(ind);
- B1(:,ind)=[];
- rho(:,:,:,ind)=[];
- pi(:,ind)=[];
-
-end
projects / valse.git / commitdiff