.RData
*.swp
*~
+/man/
+!/man/*-package.Rd
Package: valse
-Title: Variable selection with mixture of models
+Title: VAriabLe SElection with mixture of models
Date: 2016-12-01
Version: 0.1-0
Description: TODO
-Authors@R: c( person("Benjamin Auder", "Developer", role=c("ctb","cre"), email="Benjamin.Auder@math.u-psud.fr"),
- person("Benjamin Goehry", "User", role="aut", email = "Benjamin.Goehry@math.u-psud.fr"),
- person("Emilie Devijver", "User", role="ctb", email = "Emilie.Devijver@kuleuven.be"))
+Author: Benjamin Auder <Benjamin.Auder@math.u-psud.fr> [aut,cre],
+ Benjamin Goehry <Benjamin.Goehry@math.u-psud.fr> [aut]
+ Emilie Devijver <Emilie.Devijver@kuleuven.be> [aut]
+Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
Depends:
- R (>= 2.15)
-LazyData: yes
+ R (>= 3.0.0)
+Imports:
+ MASS
+Suggests:
+ parallel,
+ testthat,
+ knitr
URL: http://git.auder.net/?p=valse.git
-License: MIT
+License: MIT + file LICENSE
+VignetteBuilder: knitr
RoxygenNote: 5.0.1
--- /dev/null
+Copyright (c)
+ 2014-2016, Emilie Devijver
+ 2014-2017, Benjamin Auder
+ 2016-2017, Benjamin Goehry
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# Generated by roxygen2: do not edit by hand
-export(basic_Init_Parameters)
-export(discardSimilarModels)
-export(discardSimilarModels2)
-export(generateIO)
-export(generateIOdefault)
+export(basicInitParameters)
+export(discardSimilarModels_EMGLLF)
+export(discardSimilarModels_EMGrank)
+export(generateXY)
+export(generateXYdefault)
export(gridLambda)
-export(indicesSelection)
export(initSmallEM)
export(modelSelection)
export(selectVariables)
+++ /dev/null
-#-----------------------------------------------------------------------
-#' Initialize the parameters in a basic way (zero for the conditional mean,
-#' uniform for weights, identity for covariance matrices, and uniformly distributed forthe clustering)
-#' @param n sample size
-#' @param p number of covariates
-#' @param m size of the response
-#' @param k number of clusters
-#' @return list with phiInit, rhoInit,piInit,gamInit
-#' @export
-#-----------------------------------------------------------------------
-basic_Init_Parameters = function(n,p,m,k)
-{
- phiInit = array(0, dim=c(p,m,k))
-
- piInit = (1./k)*rep.int(1,k)
-
- rhoInit = array(0, dim=c(m,m,k))
- for(i in 1:k)
- rhoInit[,,i] = diag(m)
-
- gamInit = 0.1*array(1, dim=c(n,k))
- R = sample(1:k,n, replace=TRUE)
- for(i in 1:n)
- gamInit[i,R[i]] = 0.9
- gamInit = gamInit/sum(gamInit[1,])
-
- return (data = list(phiInit = phiInit, rhoInit = rhoInit, piInit = piInit, gamInit = gamInit))
-}
-#' Discard models which have the same relevant variables
+#' Discard models which have the same relevant variables - for EMGLLF
#'
#' @param B1 array of relevant coefficients (of size p*m*length(gridlambda))
#' @param B2 array of irrelevant coefficients (of size p*m*length(gridlambda))
#' @return a list with update B1, B2, glambda, rho and pi, and ind the vector of indices
#' of selected models.
#' @export
-discardSimilarModels = function(B1,B2,glambda,rho,pi)
+discardSimilarModels_EMGLLF = function(B1,B2,glambda,rho,pi)
{
ind = c()
for (j in 1:length(glambda))
rho = rho[,,,-ind]
pi = pi[,-ind]
- return (list(B1=B1,B2=B2,glambda=glambda,rho=rho,pi=pi,ind=ind))
+ return (list("B1"=B1,"B2"=B2,"glambda"=glambda,"rho"=rho,"pi"=pi,"ind"=ind))
+}
+
+#' Discard models which have the same relevant variables
+#' - for Lasso-rank procedure (focus on columns)
+#'
+#' @param B1 array of relevant coefficients (of size p*m*length(gridlambda))
+#' @param rho covariance matrix
+#' @param pi weight parameters
+#'
+#' @return a list with B1, in, rho, pi
+#' @export
+discardSimilarModels_EMGrank = function(B1,rho,pi)
+{
+ ind = c()
+ dim_B1 = dim(B1)
+ B2 = array(0,dim=c(dim_B1[1],dim_B1[2],dim_B1[3]))
+ sizeLambda=dim_B1[3]
+ glambda = rep(0,sizeLambda)
+
+ suppressmodel = discardSimilarModels_EMGLLF(B1,B2,glambda,rho,pi)
+ return (list("B1" = suppressmodel$B1, "ind" = suppressmodel$ind,
+ "rho" = suppressmodel$rho, "pi" = suppressmodel$pi))
}
+++ /dev/null
-#' Similar to discardSimilarModels, for Lasso-rank procedure (focus on columns)
-#'
-#' @param B1 array of relevant coefficients (of size p*m*length(gridlambda))
-#' @param rho covariance matrix
-#' @param pi weight parameters
-#'
-#' @return a list with B1, in, rho, pi
-#' @export
-#'
-discardSimilarModels2 = function(B1,rho,pi)
-{ ind = c()
- dim_B1 = dim(B1)
- B2 = array(0,dim=c(dim_B1[1],dim_B1[2],dim_B1[3]))
- sizeLambda=dim_B1[3]
- glambda = rep(0,sizeLambda)
-
- suppressmodel = discardSimilarModels(B1,B2,glambda,rho,pi)
- return (list(B1 = suppressmodel$B1, ind = suppressmodel$ind,
- rho = suppressmodel$rho, pi = suppressmodel$pi))
-}
+++ /dev/null
-#' Generate a sample of (X,Y) of size n
-#' @param covX covariance for covariates (of size p*p*K)
-#' @param covY covariance for the response vector (of size m*m*K)
-#' @param pi proportion for each cluster
-#' @param beta regression matrix
-#' @param n sample size
-#'
-#' @return list with X and Y
-#' @export
-#-----------------------------------------------------------------------
-generateIO = function(covX, covY, pi, beta, n)
-{
- p = dim(covX)[1]
-
- m = dim(covY)[1]
- k = dim(covY)[3]
-
- Y = matrix(0,n,m)
- require(mvtnorm)
- X = rmvnorm(n, mean = rep(0,p), sigma = covX)
-
- require(MASS) #simulate from a multivariate normal distribution
- for (i in 1:n)
- {
-
- for (r in 1:k)
- {
- BXir = rep(0,m)
- for (mm in 1:m)
- BXir[mm] = X[i,] %*% beta[,mm,r]
- Y[i,] = Y[i,] + pi[r] * mvrnorm(1,BXir, covY[,,r])
- }
- }
-
- return (list(X=X,Y=Y))
-}
+++ /dev/null
-#' Generate a sample of (X,Y) of size n with default values
-#' @param n sample size
-#' @param p number of covariates
-#' @param m size of the response
-#' @param k number of clusters
-#' @return list with X and Y
-#' @export
-#-----------------------------------------------------------------------
-generateIOdefault = function(n, p, m, k)
-{
- covX = diag(p)
- covY = array(0, dim=c(m,m,k))
- for(r in 1:k)
- {
- covY[,,r] = diag(m)
- }
-
- pi = rep(1./k,k)
-
- beta = array(0, dim=c(p,m,k))
- for(j in 1:p)
- {
- nonZeroCount = ceiling(m * runif(1))
- beta[j,1:nonZeroCount,] = matrix(runif(nonZeroCount*k), ncol=k)
- }
-
- sample_IO = generateIO(covX, covY, pi, beta, n)
- return (list(X=sample_IO$X,Y=sample_IO$Y))
-}
--- /dev/null
+#' Generate a sample of (X,Y) of size n
+#' @param meanX matrix of group means for covariates (of size p*K)
+#' @param covX covariance for covariates (of size p*p*K)
+#' @param covY covariance for the response vector (of size m*m*K)
+#' @param pi proportion for each cluster
+#' @param beta regression matrix
+#' @param n sample size
+#'
+#' @return list with X and Y
+#' @export
+generateXY = function(meanX, covX, covY, pi, beta, n)
+{
+ p = dim(covX)[1]
+ m = dim(covY)[1]
+ k = dim(covX)[3]
+
+ X = matrix(nrow=n,ncol=p)
+ Y = matrix(nrow=n,ncol=m)
+
+ require(MASS) #simulate from a multivariate normal distribution
+ for (i in 1:n)
+ {
+ class = sample(1:k, 1, prob=pi)
+ X[i,] = mvrnorm(1, meanX[,class], covX[,,class])
+ Y[i,] = mvrnorm(1, X[i,] %*% beta[,,class], covY[,,class])
+ }
+
+ return (list(X=X,Y=Y))
+}
+
+#' Generate a sample of (X,Y) of size n with default values
+#' @param n sample size
+#' @param p number of covariates
+#' @param m size of the response
+#' @param k number of clusters
+#' @return list with X and Y
+#' @export
+generateXYdefault = function(n, p, m, k)
+{
+ rangeX = 100
+ meanX = rangeX * matrix(1 - 2*runif(p*k), ncol=k)
+ covX = array(dim=c(p,p,k))
+ covY = array(dim=c(m,m,k))
+ for(r in 1:k)
+ {
+ covX[,,r] = diag(p)
+ covY[,,r] = diag(m)
+ }
+ pi = rep(1./k,k)
+ #initialize beta to a random number of non-zero random value
+ beta = array(0, dim=c(p,m,k))
+ for (j in 1:p)
+ {
+ nonZeroCount = sample(1:m, 1)
+ beta[j,1:nonZeroCount,] = matrix(runif(nonZeroCount*k), ncol=k)
+ }
+
+ sample_IO = generateXY(meanX, covX, covY, pi, beta, n)
+ return (list(X=sample_IO$X,Y=sample_IO$Y))
+}
+
+#' Initialize the parameters in a basic way (zero for the conditional mean, uniform for weights,
+#' identity for covariance matrices, and uniformly distributed for the clustering)
+#' @param n sample size
+#' @param p number of covariates
+#' @param m size of the response
+#' @param k number of clusters
+#' @return list with phiInit, rhoInit,piInit,gamInit
+#' @export
+basicInitParameters = function(n,p,m,k)
+{
+ phiInit = array(0, dim=c(p,m,k))
+
+ piInit = (1./k)*rep(1,k)
+
+ rhoInit = array(dim=c(m,m,k))
+ for (i in 1:k)
+ rhoInit[,,i] = diag(m)
+
+ gamInit = 0.1 * matrix(1, nrow=n, ncol=k)
+ R = sample(1:k, n, replace=TRUE)
+ for (i in 1:n)
+ gamInit[i,R[i]] = 0.9
+ gamInit = gamInit/sum(gamInit[1,])
+
+ return (list("phiInit" = phiInit, "rhoInit" = rhoInit, "piInit" = piInit, "gamInit" = gamInit))
+}
+++ /dev/null
-#' Construct the set of relevant indices -> ED: je crois que cette fonction n'est pas utile
-#'
-#' @param phi regression matrix, of size p*m
-#' @param thresh threshold to say a cofficient is equal to zero
-#'
-#' @return a list with A, a matrix with relevant indices (size = p*m) and B, a
-#' matrix with irrelevant indices (size = p*m)
-#' @export
-indicesSelection = function(phi, thresh = 1e-6)
-{
- dim_phi = dim(phi)
- p = dim_phi[1]
- m = dim_phi[2]
-
- A = matrix(0, p, m)
- B = matrix(0, p, m)
-
- for(j in 1:p)
- {
- cpt1 = 0
- cpt2 = 0
- for(mm in 1:m)
- {
- if(max(phi[j,mm,]) > thresh)
- {
- cpt1 = cpt1 + 1
- A[j,cpt] = mm
- } else
- {
- cpt2 = cpt2+1
- B[j, cpt2] = mm
- }
- }
- }
- return (list(A=A,B=B))
-}
m = ncol(Y)
p = ncol(X)
- Zinit1 = array(0, dim=c(n,20))
+ Zinit1 = array(0, dim=c(n,20))
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))
Z = Zinit1[,repet]
Z_indice = seq_len(n)[Z == r] #renvoit les indices où Z==r
- betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*% crossprod(X[Z_indice,], Y[Z_indice,])
+ betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*%
+ crossprod(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])
miniInit = 10
maxiInit = 11
- new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,tau)
+ new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],
+ gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,tau)
LLFEessai = new_EMG$LLF
LLFinit1[repet] = LLFEessai[length(LLFEessai)]
}
#' Among a collection of models, this function constructs a subcollection of models with
-#' models having strictly different dimensions, keeping the model which minimizes
+#' models having strictly different dimensions, keeping the model which minimizes
#' the likelihood if there were several with the same dimension
#'
#' @param LLF a matrix, the first column corresponds to likelihoods for several models
#' the second column corresponds to the dimensions of the corresponding models.
#'
-#' @return a list with indices, a vector of indices selected models,
+#' @return a list with indices, a vector of indices selected models,
#' and D1, a vector of corresponding dimensions
#' @export
#'
return (list(indices=indices,D1=D1))
}
+
+#TODO:
+## Programme qui sélectionne un modèle
+## proposer à l'utilisation différents critères (BIC, AIC, slope heuristic)
+++ /dev/null
-## Programme qui sélectionne un modèle
-## proposer à l'utilisation différents critères (BIC, AIC, slope heuristic)
\ No newline at end of file
#!/bin/sh
rm -f src/*.so
-rm -f adapters/*.o
-rm -f sources/*.o
+rm -f src/adapters/*.o
+rm -f src/sources/*.o
+cd src/test && make cclean && cd ../..
Trouver un jeu de données (+) intéressant (que les autres) ?
Ajouter toy datasets pour les tests (ou piocher dans MASS ?)
-ED : j'ai simulé un truc basique via 'generateIOdefault(10,5,6,2)'
+ED : j'ai simulé un truc basique via 'generateXYdefault(10,5,6,2)'
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/basicInitParameters.R
-\name{basic_Init_Parameters}
-\alias{basic_Init_Parameters}
-\title{Initialize the parameters in a basic way (zero for the conditional mean,
-uniform for weights, identity for covariance matrices, and uniformly distributed forthe clustering)}
-\usage{
-basic_Init_Parameters(n, p, m, k)
-}
-\arguments{
-\item{n}{sample size}
-
-\item{p}{number of covariates}
-
-\item{m}{size of the response}
-
-\item{k}{number of clusters}
-}
-\value{
-list with phiInit, rhoInit,piInit,gamInit
-}
-\description{
-Initialize the parameters in a basic way (zero for the conditional mean,
-uniform for weights, identity for covariance matrices, and uniformly distributed forthe clustering)
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/discardSimilarModels.R
-\name{discardSimilarModels}
-\alias{discardSimilarModels}
-\title{Discard models which have the same relevant variables}
-\usage{
-discardSimilarModels(B1, B2, glambda, rho, pi)
-}
-\arguments{
-\item{B1}{array of relevant coefficients (of size p*m*length(gridlambda))}
-
-\item{B2}{array of irrelevant coefficients (of size p*m*length(gridlambda))}
-
-\item{glambda}{grid of regularization parameters (vector)}
-
-\item{rho}{covariance matrix (of size m*m*K*size(gridLambda))}
-
-\item{pi}{weight parameters (of size K*size(gridLambda))}
-}
-\value{
-a list with update B1, B2, glambda, rho and pi, and ind the vector of indices
-of selected models.
-}
-\description{
-Discard models which have the same relevant variables
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/discardSimilarModels2.R
-\name{discardSimilarModels2}
-\alias{discardSimilarModels2}
-\title{Similar to discardSimilarModels, for Lasso-rank procedure (focus on columns)}
-\usage{
-discardSimilarModels2(B1, rho, pi)
-}
-\arguments{
-\item{B1}{array of relevant coefficients (of size p*m*length(gridlambda))}
-
-\item{rho}{covariance matrix}
-
-\item{pi}{weight parameters}
-}
-\value{
-a list with B1, in, rho, pi
-}
-\description{
-Similar to discardSimilarModels, for Lasso-rank procedure (focus on columns)
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/generateIO.R
-\name{generateIO}
-\alias{generateIO}
-\title{Generate a sample of (X,Y) of size n}
-\usage{
-generateIO(covX, covY, pi, beta, n)
-}
-\arguments{
-\item{covX}{covariance for covariates (of size p*p*K)}
-
-\item{covY}{covariance for the response vector (of size m*m*K)}
-
-\item{pi}{proportion for each cluster}
-
-\item{beta}{regression matrix}
-
-\item{n}{sample size}
-}
-\value{
-list with X and Y
-}
-\description{
-Generate a sample of (X,Y) of size n
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/generateIOdefault.R
-\name{generateIOdefault}
-\alias{generateIOdefault}
-\title{Generate a sample of (X,Y) of size n with default values}
-\usage{
-generateIOdefault(n, p, m, k)
-}
-\arguments{
-\item{n}{sample size}
-
-\item{p}{number of covariates}
-
-\item{m}{size of the response}
-
-\item{k}{number of clusters}
-}
-\value{
-list with X and Y
-}
-\description{
-Generate a sample of (X,Y) of size n with default values
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/gridLambda.R
-\name{gridLambda}
-\alias{gridLambda}
-\title{Construct the data-driven grid for the regularization parameters used for the Lasso estimator}
-\usage{
-gridLambda(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, maxi, tau)
-}
-\arguments{
-\item{phiInit}{value for phi}
-
-\item{piInit}{value for pi}
-
-\item{gamInit}{value for gamma}
-
-\item{mini}{minimum number of iterations in EM algorithm}
-
-\item{maxi}{maximum number of iterations in EM algorithm}
-
-\item{tau}{threshold to stop EM algorithm}
-
-\item{rhoInt}{value for rho}
-}
-\value{
-the grid of regularization parameters
-}
-\description{
-Construct the data-driven grid for the regularization parameters used for the Lasso estimator
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/indicesSelection.R
-\name{indicesSelection}
-\alias{indicesSelection}
-\title{Construct the set of relevant indices -> ED: je crois que cette fonction n'est pas utile}
-\usage{
-indicesSelection(phi, thresh = 1e-06)
-}
-\arguments{
-\item{phi}{regression matrix, of size p*m}
-
-\item{thresh}{threshold to say a cofficient is equal to zero}
-}
-\value{
-a list with A, a matrix with relevant indices (size = p*m) and B, a
- matrix with irrelevant indices (size = p*m)
-}
-\description{
-Construct the set of relevant indices -> ED: je crois que cette fonction n'est pas utile
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/initSmallEM.R
-\name{initSmallEM}
-\alias{initSmallEM}
-\title{initialization of the EM algorithm}
-\usage{
-initSmallEM(k, X, Y, tau)
-}
-\arguments{
-\item{k}{number of components}
-
-\item{X}{matrix of covariates (of size n*p)}
-
-\item{Y}{matrix of responses (of size n*m)}
-
-\item{tau}{threshold to stop EM algorithm}
-}
-\value{
-a list with phiInit, rhoInit, piInit, gamInit
-}
-\description{
-initialization of the EM algorithm
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/modelSelection.R
-\name{modelSelection}
-\alias{modelSelection}
-\title{Among a collection of models, this function constructs a subcollection of models with
-models having strictly different dimensions, keeping the model which minimizes
-the likelihood if there were several with the same dimension}
-\usage{
-modelSelection(LLF)
-}
-\arguments{
-\item{LLF}{a matrix, the first column corresponds to likelihoods for several models
-the second column corresponds to the dimensions of the corresponding models.}
-}
-\value{
-a list with indices, a vector of indices selected models,
- and D1, a vector of corresponding dimensions
-}
-\description{
-Among a collection of models, this function constructs a subcollection of models with
-models having strictly different dimensions, keeping the model which minimizes
-the likelihood if there were several with the same dimension
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/selectVariables.R
-\name{selectVariables}
-\alias{selectVariables}
-\title{selectVaribles
-It is a function which construct, for a given lambda, the sets of
-relevant variables and irrelevant variables.}
-\usage{
-selectVariables(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, glambda,
- X, Y, thres, tau)
-}
-\arguments{
-\item{phiInit}{an initial estimator for phi (size: p*m*k)}
-
-\item{rhoInit}{an initial estimator for rho (size: m*m*k)}
-
-\item{piInit}{an initial estimator for pi (size : k)}
-
-\item{gamInit}{an initial estimator for gamma}
-
-\item{mini}{minimum number of iterations in EM algorithm}
-
-\item{maxi}{maximum number of iterations in EM algorithm}
-
-\item{gamma}{power in the penalty}
-
-\item{glambda}{grid of regularization parameters}
-
-\item{X}{matrix of regressors}
-
-\item{Y}{matrix of responses}
-
-\item{thres}{threshold to consider a coefficient to be equal to 0}
-
-\item{tau}{threshold to say that EM algorithm has converged}
-}
-\value{
-TODO
-}
-\description{
-selectVaribles
-It is a function which construct, for a given lambda, the sets of
-relevant variables and irrelevant variables.
-}
-\examples{
-TODO
-
-}
-
+++ /dev/null
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/vec_bin.R
-\name{vec_bin}
-\alias{vec_bin}
-\title{A function needed in initSmallEM}
-\usage{
-vec_bin(X, r)
-}
-\arguments{
-\item{X}{vector with integer values}
-
-\item{r}{integer}
-}
-\value{
-a list with Z (a binary vector of size the size of X) and indices where Z is equal to 1
-}
-\description{
-A function needed in initSmallEM
-}
-
//S is already allocated, and doesn't need to be 'zeroed'
//Other local variables
- //NOTE: variables order is always [maxi],n,p,m,k
Real* gam = (Real*)malloc(n*k*sizeof(Real));
copyArray(gamInit, gam, n*k);
Real* b = (Real*)malloc(k*sizeof(Real));
Real* Gam = (Real*)malloc(n*k*sizeof(Real));
Real* X2 = (Real*)malloc(n*p*k*sizeof(Real));
Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real));
+ Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
gsl_matrix* matrix = gsl_matrix_alloc(m, m);
gsl_permutation* permutation = gsl_permutation_alloc(m);
Real* YiRhoR = (Real*)malloc(m*sizeof(Real));
Real dist = 0.;
Real dist2 = 0.;
int ite = 0;
- Real EPS = 1e-15;
- Real* dotProducts = (Real*)malloc(k*sizeof(Real));
+ const Real EPS = 1e-15;
+ const Real gaussConstM = pow(2.*M_PI,m/2.);
while (ite < mini || (ite < maxi && (dist >= tau || dist2 >= sqrt(tau))))
{
{
for (int mm=0; mm<m; mm++)
{
- //Y2(:,mm,r)=sqrt(gam(:,r)).*transpose(Y(mm,:));
+ //Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]
for (int u=0; u<n; u++)
Y2[ai(u,mm,r,n,m,k)] = sqrt(gam[mi(u,r,n,k)]) * Y[mi(u,mm,m,n)];
}
for (int i=0; i<n; i++)
{
- //X2(i,:,r)=X(i,:).*sqrt(gam(i,r));
+ //X2[i,,r] = sqrt(gam[i,r]) * X[i,]
for (int u=0; u<p; u++)
X2[ai(i,u,r,n,p,k)] = sqrt(gam[mi(i,r,n,k)]) * X[mi(i,u,n,p)];
}
for (int mm=0; mm<m; mm++)
{
- //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r);
+ //ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
for (int u=0; u<p; u++)
{
Real dotProduct = 0.;
{
for (int s=0; s<p; s++)
{
- //Gram2(j,s,r)=transpose(X2(:,j,r))*(X2(:,s,r));
+ //Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
Real dotProduct = 0.;
for (int u=0; u<n; u++)
dotProduct += X2[ai(u,j,r,n,p,k)] * X2[ai(u,s,r,n,p,k)];
// Pour pi
for (int r=0; r<k; r++)
{
- //b(r) = sum(sum(abs(phi(:,:,r))));
+ //b[r] = sum(abs(phi[,,r]))
Real sumAbsPhi = 0.;
for (int u=0; u<p; u++)
for (int v=0; v<m; v++)
sumAbsPhi += fabs(phi[ai(u,v,r,p,m,k)]);
b[r] = sumAbsPhi;
}
- //gam2 = sum(gam,1);
+ //gam2 = colSums(gam)
for (int u=0; u<k; u++)
{
Real sumOnColumn = 0.;
sumOnColumn += gam[mi(v,u,n,k)];
gam2[u] = sumOnColumn;
}
- //a=sum(gam*transpose(log(pi)));
+ //a = sum(gam %*% log(pi))
Real a = 0.;
for (int u=0; u<n; u++)
{
Real invN = 1./n;
while (!pi2AllPositive)
{
- //pi2(:)=pi(:)+0.1^kk*(1/n*gam2(:)-pi(:));
+ //pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
+ Real pow_01_kk = pow(0.1,kk);
for (int r=0; r<k; r++)
- pi2[r] = pi[r] + pow(0.1,kk) * (invN*gam2[r] - pi[r]);
+ pi2[r] = pi[r] + pow_01_kk * (invN*gam2[r] - pi[r]);
+ //pi2AllPositive = all(pi2 >= 0)
pi2AllPositive = 1;
for (int r=0; r<k; r++)
{
kk++;
}
- //t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
//(pi.^gamma)*b
Real piPowGammaDotB = 0.;
for (int v=0; v<k; v++)
Real prodGam2logPi2 = 0.;
for (int v=0; v<k; v++)
prodGam2logPi2 += gam2[v] * log(pi2[v]);
+ //t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
while (-invN*a + lambda*piPowGammaDotB < -invN*prodGam2logPi2 + lambda*pi2PowGammaDotB
&& kk<1000)
{
- //pi2=pi+0.1^kk*(1/n*gam2-pi);
+ Real pow_01_kk = pow(0.1,kk);
+ //pi2 = pi + 0.1^kk * (1/n*gam2 - pi)
for (int v=0; v<k; v++)
- pi2[v] = pi[v] + pow(0.1,kk) * (invN*gam2[v] - pi[v]);
+ pi2[v] = pi[v] + pow_01_kk * (invN*gam2[v] - pi[v]);
//pi2 was updated, so we recompute pi2PowGammaDotB and prodGam2logPi2
pi2PowGammaDotB = 0.;
for (int v=0; v<k; v++)
kk++;
}
Real t = pow(0.1,kk);
- //sum(pi+t*(pi2-pi))
+ //sum(pi + t*(pi2-pi))
Real sumPiPlusTbyDiff = 0.;
for (int v=0; v<k; v++)
sumPiPlusTbyDiff += (pi[v] + t*(pi2[v] - pi[v]));
- //pi=(pi+t*(pi2-pi))/sum(pi+t*(pi2-pi));
+ //pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi))
for (int v=0; v<k; v++)
pi[v] = (pi[v] + t*(pi2[v] - pi[v])) / sumPiPlusTbyDiff;
for (int i=0; i<n; i++)
{
//< X2(i,:,r) , phi(:,mm,r) >
- Real dotProduct = 0.0;
+ Real dotProduct = 0.;
for (int u=0; u<p; u++)
dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)];
- //ps1(i,mm,r)=Y2(i,mm,r)*dot(X2(i,:,r),phi(:,mm,r));
+ //ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
ps1[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * dotProduct;
nY21[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)];
}
- //ps(mm,r)=sum(ps1(:,mm,r));
- Real sumPs1 = 0.0;
+ //ps[mm,r] = sum(ps1[,mm,r])
+ Real sumPs1 = 0.;
for (int u=0; u<n; u++)
sumPs1 += ps1[ai(u,mm,r,n,m,k)];
ps[mi(mm,r,m,k)] = sumPs1;
- //nY2(mm,r)=sum(nY21(:,mm,r));
- Real sumNy21 = 0.0;
+ //nY2[mm,r] = sum(nY21[,mm,r])
+ Real sumNy21 = 0.;
for (int u=0; u<n; u++)
sumNy21 += nY21[ai(u,mm,r,n,m,k)];
nY2[mi(mm,r,m,k)] = sumNy21;
- //rho(mm,mm,r)=((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r)));
+ //rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
rho[ai(mm,mm,r,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
- + 4*nY2[mi(mm,r,m,k)] * (gam2[r]) ) ) / (2*nY2[mi(mm,r,m,k)]);
+ + 4*nY2[mi(mm,r,m,k)] * gam2[r] ) ) / (2*nY2[mi(mm,r,m,k)]);
}
}
for (int r=0; r<k; r++)
{
for (int mm=0; mm<m; mm++)
{
- //sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))+sum(phi(j+1:p,mm,r)
- // .*transpose(Gram2(j,j+1:p,r)))
+ //sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r])
Real dotPhiGram2 = 0.0;
for (int u=0; u<j; u++)
dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
+ //sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r])
for (int u=j+1; u<p; u++)
dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
- //S(j,r,mm)=-rho(mm,mm,r)*ps2(j,mm,r)+sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))
- // +sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r)));
+ //S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] +
+ // (if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
+ // (if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] + dotPhiGram2;
- if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow(pi[r],gamma))
+ Real pow_pir_gamma = pow(pi[r],gamma);
+ if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow_pir_gamma)
phi[ai(j,mm,r,p,m,k)] = 0;
- else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow(pi[r],gamma))
- phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow(pi[r],gamma) - S[ai(j,mm,r,p,m,k)])
+ else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow_pir_gamma)
+ {
+ phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow_pir_gamma - S[ai(j,mm,r,p,m,k)])
/ Gram2[ai(j,j,r,p,p,k)];
+ }
else
- phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow(pi[r],gamma) + S[ai(j,mm,r,p,m,k)])
+ {
+ phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow_pir_gamma + S[ai(j,mm,r,p,m,k)])
/ Gram2[ai(j,j,r,p,p,k)];
+ }
}
}
}
Real sumLogLLF2 = 0.0;
for (int i=0; i<n; i++)
{
- Real sumLLF1 = 0.0;
- Real sumGamI = 0.0;
- Real minDotProduct = INFINITY;
+ Real minSqNorm2 = INFINITY;
for (int r=0; r<k; r++)
{
- //Compute
- //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
- // *phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
- //split in several sub-steps
-
- //compute Y(i,:)*rho(:,:,r)
+ //compute Y[i,]%*%rho[,,r]
for (int u=0; u<m; u++)
{
YiRhoR[u] = 0.0;
XiPhiR[u] += X[mi(i,v,n,p)] * phi[ai(v,u,r,p,m,k)];
}
- //compute dotProduct
- // < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
- dotProducts[r] = 0.0;
+ //compute sq norm || Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) ||_2^2
+ sqNorm2[r] = 0.0;
for (int u=0; u<m; u++)
- dotProducts[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
- if (dotProducts[r] < minDotProduct)
- minDotProduct = dotProducts[r];
+ sqNorm2[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
+ if (sqNorm2[r] < minSqNorm2)
+ minSqNorm2 = sqNorm2[r];
}
- Real shift = 0.5*minDotProduct;
+ Real shift = 0.5*minSqNorm2;
+
+ Real sumLLF1 = 0.0;
+ Real sumGamI = 0.0;
for (int r=0; r<k; r++)
{
- //compute det(rho(:,:,r)) [TODO: avoid re-computations]
+ //compute det(rho[,,r]) [TODO: avoid re-computations]
for (int u=0; u<m; u++)
{
for (int v=0; v<m; v++)
gsl_linalg_LU_decomp(matrix, permutation, &signum);
Real detRhoR = gsl_linalg_LU_det(matrix, signum);
- Gam[mi(i,r,n,k)] = pi[r] * detRhoR * exp(-0.5*dotProducts[r] + shift);
- sumLLF1 += Gam[mi(i,r,n,k)] / pow(2*M_PI,m/2.0);
+ //FIXME: det(rho[,,r]) too small(?!). See EMGLLF.R
+ Gam[mi(i,r,n,k)] = pi[r] * exp(-0.5*sqNorm2[r] + shift) * detRhoR;
+ sumLLF1 += Gam[mi(i,r,n,k)] / gaussConstM;
sumGamI += Gam[mi(i,r,n,k)];
}
sumLogLLF2 += log(sumLLF1);
for (int r=0; r<k; r++)
{
- //gam(i,r)=Gam(i,r)/sum(Gam(i,:));
- gam[mi(i,r,n,k)] = sumGamI > EPS
- ? Gam[mi(i,r,n,k)] / sumGamI
- : 0.0;
+ //gam[i,] = Gam[i,] / sumGamI
+ gam[mi(i,r,n,k)] = sumGamI > EPS ? Gam[mi(i,r,n,k)] / sumGamI : 0.;
}
}
-
- //sum(pen(ite,:))
+
+ //sumPen = sum(pi^gamma * b)
Real sumPen = 0.0;
for (int r=0; r<k; r++)
sumPen += pow(pi[r],gamma) * b[r];
- //LLF(ite)=-1/n*sum(log(LLF2(ite,:)))+lambda*sum(pen(ite,:));
+ //LLF[ite] = -sumLogLLF2/n + lambda*sumPen
LLF[ite] = -invN * sumLogLLF2 + lambda * sumPen;
- if (ite == 0)
- dist = LLF[ite];
- else
- dist = (LLF[ite] - LLF[ite-1]) / (1.0 + fabs(LLF[ite]));
-
- //Dist1=max(max((abs(phi-Phi))./(1+abs(phi))));
+ dist = ite==0 ? LLF[ite] : (LLF[ite] - LLF[ite-1]) / (1.0 + fabs(LLF[ite]));
+
+ //Dist1 = max( abs(phi-Phi) / (1+abs(phi)) )
Real Dist1 = 0.0;
for (int u=0; u<p; u++)
{
}
}
}
- //Dist2=max(max((abs(rho-Rho))./(1+abs(rho))));
+ //Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) )
Real Dist2 = 0.0;
for (int u=0; u<m; u++)
{
}
}
}
- //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi))));
+ //Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)))
Real Dist3 = 0.0;
for (int u=0; u<n; u++)
{
dist2 = Dist2;
if (Dist3 > dist2)
dist2 = Dist3;
-
+
ite++;
}
-
+
//free memory
free(b);
free(gam);
free(pi2);
free(X2);
free(Y2);
- free(dotProducts);
+ free(sqNorm2);
}
// TODO: comment on constructionModelesLassoRank purpose
void constructionModelesLassoRank_core(
// IN parameters
- const Real* Pi,// parametre initial des proportions
- const Real* Rho, // parametre initial de variance renormalisé
+ const Real* pi,// parametre initial des proportions
+ const Real* rho, // parametre initial de variance renormalisé
int mini, // nombre minimal d'itérations dans l'algorithme EM
int maxi, // nombre maximal d'itérations dans l'algorithme EM
const Real* X,// régresseurs
int deltaRank = rangmax-rangmin+1;
int Size = (int)pow(deltaRank,k);
int* Rank = (int*)malloc(Size*k*sizeof(int));
-for (int r=0; r<k; r++)
-{
- //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.
+ for (int r=0; r<k; r++)
+ {
+ // 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 :
+ // ça remplit la colonne
+ // Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre,
+ // et on fait ça (rangmax-rangmin)^2 fois
+ // ...
+ // Dans la dernière, on répète chaque chiffre une fois,
+ // et on fait ça (rangmin-rangmax)^(k-1) fois.
int indexInRank = 0;
int value = 0;
while (indexInRank < Size)
{
- for (int u=0; u<pow(deltaRank,k-r-1); u++)
+ int upperBoundU = (int)pow(deltaRank,k-r-1);
+ for (int u=0; u<upperBoundU; u++)
Rank[mi(indexInRank++,r,Size,k)] = rangmin + value;
value = (value+1) % deltaRank;
}
if (A1[mi(j,lambdaIndex,p,L)] != 0)
active[longueurActive++] = A1[mi(j,lambdaIndex,p,L)] - 1;
}
-
if (longueurActive == 0)
continue;
Real LLF;
for (int j=0; j<Size; j++)
{
- //[phiLambda,LLF] = EMGrank(Pi(:,lambdaIndex),Rho(:,:,:,lambdaIndex),mini,maxi,X(:,active),Y,tau,Rank(j,:));
int* rank = (int*)malloc(k*sizeof(int));
for (int r=0; r<k; r++)
rank[r] = Rank[mi(j,r,Size,k)];
for (int jj=0; jj<longueurActive; jj++)
Xactive[mi(i,jj,n,longueurActive)] = X[mi(i,active[jj],n,p)];
}
- Real* PiLambda = (Real*)malloc(k*sizeof(Real));
+ Real* piLambda = (Real*)malloc(k*sizeof(Real));
for (int r=0; r<k; r++)
- PiLambda[r] = Pi[mi(r,lambdaIndex,k,L)];
- Real* RhoLambda = (Real*)malloc(m*m*k*sizeof(Real));
+ piLambda[r] = pi[mi(r,lambdaIndex,k,L)];
+ Real* rhoLambda = (Real*)malloc(m*m*k*sizeof(Real));
for (int u=0; u<m; u++)
{
for (int v=0; v<m; v++)
{
for (int r=0; r<k; r++)
- RhoLambda[ai(u,v,r,m,m,k)] = Rho[ai4(u,v,r,lambdaIndex,m,m,k,L)];
+ rhoLambda[ai(u,v,r,m,m,k)] = rho[ai4(u,v,r,lambdaIndex,m,m,k,L)];
}
}
- EMGrank_core(PiLambda,RhoLambda,mini,maxi,Xactive,Y,tau,rank,
+ EMGrank_core(piLambda,rhoLambda,mini,maxi,Xactive,Y,tau,rank,
phiLambda,&LLF,
n,longueurActive,m,k);
free(rank);
free(Xactive);
- free(PiLambda);
- free(RhoLambda);
- //llh((lambdaIndex-1)*Size+j,:) = [LLF, dot(Rank(j,:), length(active)-Rank(j,:)+m)];
+ free(piLambda);
+ free(rhoLambda);
+ //llh[(lambdaIndex-1)*Size+j,] = c(LLF, ...)
llh[mi(lambdaIndex*Size+j,0,L*Size,2)] = LLF;
- //dot(Rank(j,:), length(active)-Rank(j,:)+m)
+ //sum(Rank[j,] * (length(active)- Rank[j,] + m)) )
Real dotProduct = 0.0;
for (int r=0; r<k; r++)
dotProduct += Rank[mi(j,r,Size,k)] * (longueurActive-Rank[mi(j,r,Size,k)]+m);
llh[mi(lambdaIndex*Size+j,1,Size*L,2)] = dotProduct;
- //phi(active,:,:,(lambdaIndex-1)*Size+j) = phiLambda;
+ //phi[active,,,(lambdaIndex-1)*Size+j] = res$phi
for (int jj=0; jj<longueurActive; jj++)
{
for (int mm=0; mm<m; mm++)
{
for (int r=0; r<k; r++)
- phi[ai5(active[jj],mm,r,lambdaIndex,j,p,m,k,L,Size)] = phiLambda[jj*m*k+mm*k+r];
+ {
+ phi[ai4(active[jj],mm,r,lambdaIndex*Size+j,p,m,k,L*Size)] =
+ phiLambda[ai(jj,mm,r,longueurActive,m,k)];
+ }
}
}
}
LIB_SRC = $(wildcard ../sources/*.c)
LIB_OBJ = $(LIB_SRC:.c=.o)
INCLUDES = -I../sources
+TESTS = test.EMGLLF test.EMGrank test.constructionModelesLassoMLE test.EMGrank\
+ test.constructionModelesLassoRank test.selectionTotale
-all: $(LIB) test.EMGLLF test.EMGrank test.constructionModelesLassoMLE test.EMGrank test.constructionModelesLassoRank test.selectionTotale
+all: $(LIB) $(TESTS)
$(LIB): $(LIB_OBJ)
$(CC) -shared -o $@ $^ $(LDFLAGS)
-test.EMGLLF: test.EMGLLF.o test_utils.o
+test.EMGLLF: $(LIB) test.EMGLLF.o test_utils.o
$(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS)
-test.constructionModelesLassoMLE: test.constructionModelesLassoMLE.o test_utils.o
+test.EMGrank: $(LIB) test.EMGrank.o test_utils.o
$(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS)
-test.EMGrank: test.EMGrank.o test_utils.o
+test.constructionModelesLassoMLE: $(LIB) test.constructionModelesLassoMLE.o test_utils.o
$(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS)
-test.constructionModelesLassoRank: test.constructionModelesLassoRank.o test_utils.o
+test.constructionModelesLassoRank: $(LIB) test.constructionModelesLassoRank.o test_utils.o
$(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS)
-test.selectionTotale: test.selectionTotale.o test_utils.o
+test.selectionTotale: $(LIB) test.selectionTotale.o test_utils.o
$(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS)
%.o: %.c
rm -f *.o ../sources/*.o
cclean: clean
- rm -f $(LIB) $(TEST)
+ rm -f $(LIB) $(TESTS)
.PHONY: all clean cclean
dir.create(testFolder, showWarnings=FALSE, mode="0755")
require(valse)
- params = valse:::basic_Init_Parameters(n, p, m, k)
- io = generateIOdefault(n, p, m, k)
+ params = valse:::basicInitParameters(n, p, m, k)
+ xy = valse:::generateXYdefault(n, p, m, k)
#save inputs
write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
row.names=F, col.names=F)
write.table(as.double(lambda), paste(testFolder,"lambda",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+ write.table(as.double(xy$X), paste(testFolder,"X",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+ write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
row.names=F, col.names=F)
write.table(as.double(tau), paste(testFolder,"tau",sep=""),
row.names=F, col.names=F)
- write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+ write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensxyns",sep=""),
row.names=F, col.names=F)
res = EMGLLF(params$phiInit,params$rhoInit,params$piInit,params$gamInit,mini,maxi,
- gamma,lambda,io$X,io$Y,tau)
+ gamma,lambda,xy$X,xy$Y,tau)
#save outputs
write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F, col.names=F)
generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0,
rank = c(1,2,4))
{
- testFolder = "data/"
- dir.create(testFolder, showWarnings=FALSE, mode="0755")
-
tau = 1e-6
pi = rep(1.0/k, k)
- rho = array(0, dim=c(m,m,k))
- for(i in 1:k){
+ rho = array(dim=c(m,m,k))
+ for(i in 1:k)
rho[,,i] = diag(1,m)
- }
require(valse)
- io = valse:::generateIOdefault(n, p, m, k)
+ xy = valse:::generateXYdefault(n, p, m, k)
+ testFolder = "data/"
+ dir.create(testFolder, showWarnings=FALSE, mode="0755")
#save inputs
write.table(as.double(rho), paste(testFolder,"rho",sep=""),
row.names=F, col.names=F)
row.names=F, col.names=F)
write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+ write.table(as.double(xy$X), paste(testFolder,"X",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+ write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
row.names=F, col.names=F)
write.table(as.double(tau), paste(testFolder,"tau",sep=""),
row.names=F, col.names=F)
write.table(as.integer(rank), paste(testFolder,"rank",sep=""),
row.names=F, col.names=F)
- write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+ write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensxyns",sep=""),
row.names=F, col.names=F)
- res = EMGrank(pi,rho,mini,maxi,io$X,io$Y,tau,rank)
+ res = EMGrank(pi,rho,mini,maxi,xy$X,xy$Y,tau,rank)
#save output
write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F,col.names=F)
+++ /dev/null
-function[] = generateRunSaveTest_EMGrank(n, p, m, k, mini, maxi, gamma, rank, varargin)
-
- %set defaults for optional inputs
- optargs = {200 15 10 3 5 10 1.0 1:3};
- %replace defaults by user parameters
- optargs(1:length(varargin)) = varargin;
- [n, p, m, k, mini, maxi, gamma, rank] = optargs{:};
- mini = int64(mini);
- maxi = int64(maxi);
- rank = int64(rank);
- tau = 1e-6;
-
- Pi = (1.0/k)*ones(1,k);
- Rho = zeros(m,m,k);
- for r=1:k
- Rho(:,:,r) = eye(m);
- end
-
- %Generate X and Y
- [X, Y, ~] = generateIOdefault(n, p, m, k);
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- testFolder = 'data/';
- mkdir(testFolder);
- delimiter = ' ';
-
- %save inputs
- dlmwrite(strcat(testFolder,'Rho'), reshape(Rho,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Pi'), Pi, delimiter);
- dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
- dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
- dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
- dlmwrite(strcat(testFolder,'rank'), rank, delimiter);
- dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k], delimiter);;
-
- [phi,LLF] = EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank);
-
- %save output
- dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'LLF'), LLF, delimiter);
-
-end
generateRunSaveTest_constructionModelesLassoMLE = function(n=200, p=15, m=10, k=3, mini=5,
maxi=10, gamma=1.0, glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99))
{
- testFolder = "data/"
- dir.create(testFolder, showWarnings=FALSE, mode="0755")
-
- require(valse)
- params = valse:::basic_Init_Parameters(n, p, m, k)
- A2 = array(0, dim=c(p, m+1, L))
+ tau = 1e-6;
+ seuil = 1e-15;
+ L = length(glambda);
A1 = array(0, dim=c(p, m+1, L))
+ A2 = array(0, dim=c(p, m+1, L))
for (i in 1:L)
{
- A2[,1,i] = seq_len(p)
- A2[,2,i] = seq_len(5)
A1[,1,i] = seq_len(p)
- A1[,2,i] = seq_len(5)
+ A1[1:5,2,i] = seq_len(5)
+ A2[,1,i] = seq_len(p)
+ A2[1:5,2,i] = seq_len(5)
}
- io = generateIOdefault(n, p, m, k)
+ require(valse)
+ params = valse:::basicInitParameters(n, p, m, k)
+ xy = valse:::generateXYdefault(n, p, m, k)
+ testFolder = "data/"
+ dir.create(testFolder, showWarnings=FALSE, mode="0755")
#save inputs
write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
row.names=F, col.names=F)
row.names=F, col.names=F)
write.table(as.double(lambda), paste(testFolder,"lambda",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+ write.table(as.double(xy$X), paste(testFolder,"X",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+ write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
+ row.names=F, col.names=F)
+ write.table(as.double(seuil), paste(testFolder,"seuil",sep=""),
row.names=F, col.names=F)
write.table(as.double(tau), paste(testFolder,"tau",sep=""),
row.names=F, col.names=F)
- write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+ write.table(as.double(A1), paste(testFolder,"A1",sep=""),
+ row.names=F, col.names=F)
+ write.table(as.double(A2), paste(testFolder,"A2",sep=""),
+ row.names=F, col.names=F)
+ write.table(as.integer(c(n,p,m,k,L)), paste(testFolder,"dimensions",sep=""),
row.names=F, col.names=F)
res = constructionModelesLassoMLE(
params$phiInit,params$rhoInit,params$piInit,params$gamInit,
- mini,maxi,gamma,glambda,io$X,io$Y,seuil,tau,A1,A2)
+ mini,maxi,gamma,glambda,xy$X,xy$Y,seuil,tau,A1,A2)
#save output
write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F, col.names=F)
+++ /dev/null
-function[] = generateRunSaveTest_constructionModelesLassoMLE(n, p, m, k, mini, maxi, gamma, glambda, varargin)
-
- %set defaults for optional inputs
- optargs = {200 15 10 3 5 10 1.0 [0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99]};
- %replace defaults by user parameters
- optargs(1:length(varargin)) = varargin;
- [n, p, m, k, mini, maxi, gamma, glambda] = optargs{:};
- tau = 1e-6;
- seuil = 1e-15;
- mini = int64(mini);
- maxi = int64(maxi);
- L = length(glambda);
-
- %Generate phiInit,piInit,...
- [phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n, p, m, k);
-
- %Generate X and Y
- [X, Y, ~] = generateIOdefault(n, p, m, k);
-
- A2 = zeros(p,m+1,L,'int64');
- for i=1:L
- A2(:,1,i) = 1:p;
- A2(1:5,2,i) = 1:5;
- end
- A1 = zeros(p,m+1,L,'int64');
- for i=1:L
- A1(:,1,i) = 1:p;
- A1(1:5,2,i) = 1:5;
- end
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- testFolder = 'data/';
- mkdir(testFolder);
- delimiter = ' ';
-
- %save inputs
- dlmwrite(strcat(testFolder,'phiInit'), reshape(phiInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'rhoInit'), reshape(rhoInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'piInit'), piInit, delimiter);
- dlmwrite(strcat(testFolder,'gamInit'), reshape(gamInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
- dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
- dlmwrite(strcat(testFolder,'gamma'), gamma, delimiter);
- dlmwrite(strcat(testFolder,'glambda'), glambda, delimiter);
- dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'seuil'), seuil, delimiter);
- dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
- dlmwrite(strcat(testFolder,'A1'), reshape(A1,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'A2'), reshape(A2,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
-
- [phi,rho,pi,llh] = constructionModelesLassoMLE(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2);
-
- %save output
- dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'rho'), reshape(rho,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'pi'), reshape(pi,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'llh'), reshape(llh,1,[]), delimiter);
-
-end
-generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, L=12, mini=5,
+generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, k=3, L=12, mini=5,
maxi=10, gamma=1.0, rangmin=3, rangmax=6)
{
- testFolder = "data/"
- dir.create(testFolder, showWarnings=FALSE, mode="0755")
-
tau = 1e-6
- pi = matrix(0, k,L)
- for(i in 1:L){
- pi[,i] = rep(1.0/k, k)
- }
- rho = array(0, dim=c(m,m,k,L))
- for(l in 1:L){
- for(r in 1:k){
+ pi = matrix(1./k, nrow=k, ncol=L)
+ rho = array(dim=c(m,m,k,L))
+ for (l in 1:L)
+ {
+ for (r in 1:k)
rho[,,r,l] = diag(1,m)
- }
}
+ A1 = matrix(seq_len(p), nrow=p, ncol=L)
require(valse)
- io = valse:::generateIOdefault(n, p, m, k)
- A1 = matrix(0,p,L)
- for(i in 1:L){
- A1[,i] = seq(1,p)
- }
+ xy = valse:::generateXYdefault(n, p, m, k)
+ testFolder = "data/"
+ dir.create(testFolder, showWarnings=FALSE, mode="0755")
#save inputs
- write.table(as.double(rho),paste(testFolder,"rho",sep=""),
- row.names=F, col.names=F)
write.table(as.double(pi), paste(testFolder,"pi",sep=""),
row.names=F, col.names=F)
+ write.table(as.double(rho),paste(testFolder,"rho",sep=""),
+ row.names=F, col.names=F)
write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
row.names=F, col.names=F)
write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+ write.table(as.double(xy$X), paste(testFolder,"X",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+ write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
row.names=F, col.names=F)
write.table(as.double(tau),paste(testFolder,"tau",sep=""),
row.names=F, col.names=F)
row.names=F, col.names=F)
write.table(as.integer(rangmax),paste(testFolder,"rangmax",sep=""),
row.names=F, col.names=F)
- write.table(as.integer(c(n,p,m,k)),paste(testFolder,"dimensions",sep=""),
+ write.table(as.integer(c(n,p,m,k,L)),paste(testFolder,"dimensions",sep=""),
row.names=F, col.names=F)
- res = constructionModelesLassoRank(pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
+ res = constructionModelesLassoRank(pi,rho,mini,maxi,xy$X,xy$Y,tau,A1,rangmin,rangmax)
#save output
write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F, col.names=F)
+++ /dev/null
-function[] = generateRunSaveTest_constructionModelesLassoRank(n, p, m, k, L, mini, maxi, gamma, rangmin, rangmax, varargin)
-
- %set defaults for optional inputs
- optargs = {200 15 10 3 12 5 10 1.0 3 6};
- %replace defaults by user parameters
- optargs(1:length(varargin)) = varargin;
- [n, p, m, k, L, mini, maxi, gamma, rangmin, rangmax] = optargs{:};
- mini = int64(mini);
- maxi = int64(maxi);
- rangmin = int64(rangmin);
- rangmax = int64(rangmax);
- tau = 1e-6;
-
- Pi = zeros(k,L);
- for l=1:L
- Pi(:,l) = (1.0/k)*ones(1,k);
- end
- Rho = zeros(m,m,k,L);
- for l=1:L
- for r=1:k
- Rho(:,:,r,l) = eye(m);
- end
- end
-
- %Generate X and Y
- [X, Y, ~] = generateIOdefault(n, p, m, k);
-
- A1 = zeros(p,L,'int64');
- for i=1:L
- A1(:,i) = 1:p;
- end
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- testFolder = 'data/';
- mkdir(testFolder);
- delimiter = ' ';
-
- %save inputs
- dlmwrite(strcat(testFolder,'Rho'), reshape(Rho,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Pi'), reshape(Pi,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
- dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
- dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
- dlmwrite(strcat(testFolder,'A1'), reshape(A1,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'rangmin'), rangmin, delimiter);
- dlmwrite(strcat(testFolder,'rangmax'), rangmax, delimiter);
- dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
-
- [phi,llh] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax);
-
- %save output
- Size = (rangmax-rangmin+1)^k;
- dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'llh'), reshape(llh,1,[]), delimiter);
-
-end
generateRunSaveTest_selectiontotale= function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0,
glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99))
{
- testFolder = "data/"
- dir.create(testFolder, showWarnings=FALSE, mode="0755")
-
+ tau = 1e-6;
+ seuil = 1e-15;
require(valse)
- params = valse:::basic_Init_Parameters(n, p, m, k)
- io = valse:::generateIOdefault(n, p, m, k)
+ params = valse:::basicInitParameters(n, p, m, k)
+ xy = valse:::generateXYdefault(n, p, m, k)
+ L = length(glambda)
+ testFolder = "data/"
+ dir.create(testFolder, showWarnings=FALSE, mode="0755")
#save inputs
write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
row.names=F, col.names=F)
row.names=F, col.names=F)
write.table(as.double(glambda), paste(testFolder,"lambda",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+ write.table(as.double(xy$X), paste(testFolder,"X",sep=""),
+ row.names=F, col.names=F)
+ write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
row.names=F, col.names=F)
- write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+ write.table(as.double(seuil), paste(testFolder,"seuil",sep=""),
row.names=F, col.names=F)
write.table(as.double(tau), paste(testFolder,"tau",sep=""),
row.names=F, col.names=F)
- write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+ write.table(as.integer(c(n,p,m,k,L)), paste(testFolder,"dimensions",sep=""),
row.names=F, col.names=F)
res = selectiontotale(params$phiInit,params$rhoInit,params$piInit,params$gamInit,
- mini,maxi,gamma,glambda,io$X,io$Y,seuil, tau)
+ mini,maxi,gamma,glambda,xy$X,xy$Y,seuil, tau)
#save output
write.table(as.double(res$A1), paste(testFolder,"A1",sep=""), row.names=F, col.names=F)
+++ /dev/null
-function[] = generateRunSaveTest_selectiontotale(n, p, m, k, mini, maxi, gamma, glambda, varargin)
-
- %set defaults for optional inputs
- optargs = {200 15 10 3 5 10 1.0 [0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99]};
- %replace defaults by user parameters
- optargs(1:length(varargin)) = varargin;
- [n, p, m, k, mini, maxi, gamma, glambda] = optargs{:};
- tau = 1e-6;
- seuil = 1e-15;
- mini = int64(mini);
- maxi = int64(maxi);
- L = length(glambda);
-
- %Generate phiInit,piInit,...
- [phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n, p, m, k);
-
- %Generate X and Y
- [X, Y, ~] = generateIOdefault(n, p, m, k);
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- testFolder = 'data/';
- mkdir(testFolder);
- delimiter = ' ';
-
- %save inputs
- dlmwrite(strcat(testFolder,'phiInit'), reshape(phiInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'rhoInit'), reshape(rhoInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'piInit'), piInit, delimiter);
- dlmwrite(strcat(testFolder,'gamInit'), reshape(gamInit,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
- dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
- dlmwrite(strcat(testFolder,'gamma'), gamma, delimiter);
- dlmwrite(strcat(testFolder,'glambda'), glambda, delimiter);
- dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'seuil'), seuil, delimiter);
- dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
- dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
-
- [A1,A2,Rho,Pi] = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau);
-
- %save output
- dlmwrite(strcat(testFolder,'A1'), reshape(A1,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'A2'), reshape(A2,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Rho'), reshape(Rho,1,[]), delimiter);
- dlmwrite(strcat(testFolder,'Pi'), reshape(Pi,1,[]), delimiter);
-
-end
-EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){
+EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
+{
#matrix dimensions
n = dim(X)[1]
p = dim(phiInit)[1]
#init outputs
phi = phiInit
rho = rhoInit
- Pi = piInit
+ pi = piInit
LLF = rep(0, maxi)
S = array(0, dim=c(p,m,k))
-
gam = gamInit
Gram2 = array(0, dim=c(p,p,k))
ps2 = array(0, dim=c(p,m,k))
dist = 0
dist2 = 0
ite = 1
- Pi2 = rep(0, k)
+ pi2 = rep(0, k)
ps = matrix(0, m,k)
nY2 = matrix(0, m,k)
ps1 = array(0, dim=c(n,m,k))
Gam = matrix(0, n,k)
EPS = 1E-15
- while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))){
+ while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau))))
+ {
Phi = phi
Rho = rho
- PI = Pi
+ Pi = pi
+
#calcul associé à Y et X
- for(r in 1:k){
- for(mm in 1:m){
- Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm] ##bon calcul ? idem pour X2 ??...
- }
- for(i in 1:n){
- X2[i,,r] = X[i,] *sqrt(gam[i,r])
- }
- for(mm in 1:m){
+ for(r in 1:k)
+ {
+ for (mm in 1:m)
+ Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]
+ for (i in 1:n)
+ X2[i,,r] = sqrt(gam[i,r]) * X[i,]
+ for (mm in 1:m)
ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
- }
- for(j in 1:p){
- for(s in 1:p){
+ for (j in 1:p)
+ {
+ for (s in 1:p)
Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
- }
}
}
##########
#pour pi
- for(r in 1:k){
- b[r] = sum(sum(abs(phi[,,r])))
- }
+ for (r in 1:k)
+ b[r] = sum(abs(phi[,,r]))
gam2 = colSums(gam)
- a = sum(gam%*%(log(Pi)))
+ a = sum(gam %*% log(pi))
#tant que les props sont negatives
kk = 0
pi2AllPositive = FALSE
- while(pi2AllPositive == FALSE){
- Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi)
- pi2AllPositive = TRUE
- for(r in 1:k){
- if(Pi2[r] < 0){
- pi2AllPositive = false;
- break
- }
- }
+ while (!pi2AllPositive)
+ {
+ pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi)
+ pi2AllPositive = all(pi2 >= 0)
kk = kk+1
}
- #t[m]la plus grande valeur dans la grille O.1^k tel que ce soit
- #décroissante ou constante
- while((-1/n*a+lambda*((Pi^gamma)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(b)) && kk<1000){
- Pi2 = Pi+0.1^kk*(1/n*gam2-Pi)
- kk = kk+1
+ #t[m] la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
+ while( kk < 1000 && -a/n + lambda * sum(pi^gamma * b) <
+ -sum(gam2 * log(pi2))/n + lambda * sum(pi2^gamma * b) )
+ {
+ pi2 = pi + 0.1^kk * (1/n*gam2 - pi)
+ kk = kk + 1
}
- t = 0.1^(kk)
- Pi = (Pi+t*(Pi2-Pi)) / sum(Pi+t*(Pi2-Pi))
+ t = 0.1^kk
+ pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi))
#Pour phi et rho
- for(r in 1:k){
- for(mm in 1:m){
- for(i in 1:n){
- ps1[i,mm,r] = Y2[i,mm,r] * (X2[i,,r]%*%(phi[,mm,r]))
- nY21[i,mm,r] = (Y2[i,mm,r])^2
+ for (r in 1:k)
+ {
+ for (mm in 1:m)
+ {
+ for (i in 1:n)
+ {
+ ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
+ nY21[i,mm,r] = Y2[i,mm,r]^2
}
ps[mm,r] = sum(ps1[,mm,r])
nY2[mm,r] = sum(nY21[,mm,r])
- rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r]))
+ rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r])
}
}
- for(r in 1:k){
- p1 = p-1
- for(j in 1:p1){
- for(mm in 1:m){
- j1 = j-1
- j2 = j+1
- v1 = c(1:j1)
- v2 = c(j2:p)
- S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r]) + phi[v2,mm,r]%*%(Gram2[j,v2,r]) #erreur indice
- if(abs(S[j,mm,r]) <= n*lambda*(Pi[r]^gamma)){
+ for (r in 1:k)
+ {
+ for (j in 1:p)
+ {
+ for (mm in 1:m)
+ {
+ S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] +
+ (if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
+ (if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
+ if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
phi[j,mm,r]=0
- }else{
- if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)){
- phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r]
- }else{
- phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
- }
- }
+ else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
+ phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
+ else
+ phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r]
}
}
}
-
+
##########
#Etape E #
##########
sumLogLLF2 = 0
- for(i in 1:n){
- #precompute dot products to numerically adjust their values
- dotProducts = rep(0,k)
- for(r in 1:k){
- dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])
- }
- shift = 0.5*min(dotProducts)
-
+ for (i in 1:n)
+ {
+ #precompute sq norms to numerically adjust their values
+ sqNorm2 = rep(0,k)
+ for (r in 1:k)
+ sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 )
+ shift = 0.5*min(sqNorm2)
+
#compute Gam(:,:) using shift determined above
sumLLF1 = 0.0;
- for(r in 1:k){
- Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift)
- sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2)
+ for (r in 1:k)
+ {
+ #FIXME: numerical problems, because 0 < det(Rho[,,r] < EPS; what to do ?!
+ # consequence: error in while() at line 77
+ Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r] + shift) * det(rho[,,r])
+ sumLLF1 = sumLLF1 + Gam[i,r] / (2*base::pi)^(m/2)
}
sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
sumGamI = sum(Gam[i,])
if(sumGamI > EPS)
gam[i,] = Gam[i,] / sumGamI
else
- gam[i,] = rep(0,k)
- }
-
-
- sumPen = 0
- for(r in 1:k){
- sumPen = sumPen + Pi[r]^gamma^b[r]
+ gam[i,] = rep(0,k)
}
- LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
-
- if(ite == 1)
- dist = LLF[ite]
- else
- dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
-
- Dist1=max(max(max((abs(phi-Phi))/(1+abs(phi)))))
- Dist2=max(max(max((abs(rho-Rho))/(1+abs(rho)))))
- Dist3=max(max((abs(Pi-PI))/(1+abs(PI))))
- dist2=max(c(Dist1,Dist2,Dist3))
-
- ite=ite+1
+
+ sumPen = sum(pi^gamma * b)
+ LLF[ite] = -sumLogLLF2/n + lambda*sumPen
+
+ dist = ifelse( ite == 1, LLF[ite], (LLF[ite]-LLF[ite-1]) / (1+abs(LLF[ite])) )
+
+ Dist1 = max( (abs(phi-Phi)) / (1+abs(phi)) )
+ Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) )
+ Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)) )
+ dist2 = max(Dist1,Dist2,Dist3)
+
+ ite = ite+1
}
-
- Pi = t(Pi)
- return(list("phi"=phi, "rho"=rho, "pi"=Pi, "LLF"=LLF, "S"=S))
+
+ return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S))
}
+++ /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
}
require(MASS)
-EMGrank = function(Pi, Rho, mini, maxi, X, Y, tau, rank){
+EMGrank = function(Pi, Rho, mini, maxi, X, Y, tau, rank)
+{
#matrix dimensions
n = dim(X)[1]
p = dim(X)[2]
#init outputs
phi = array(0, dim=c(p,m,k))
Z = rep(1, n)
-# Pi = piInit
LLF = 0
#local variables
Phi = array(0, dim=c(p,m,k))
- deltaPhi = c(0)
- sumDeltaPhi = 0
+ deltaPhi = c()
+ sumDeltaPhi = 0.
deltaPhiBufferSize = 20
#main loop
ite = 1
- while(ite<=mini || (ite<=maxi && sumDeltaPhi>tau))
+ while (ite<=mini || (ite<=maxi && sumDeltaPhi>tau))
{
#M step: Mise à jour de Beta (et donc phi)
for(r in 1:k)
s = svd( ginv(crossprod(matricize(X[Z_indice,]))) %*%
crossprod(matricize(X[Z_indice,]),matricize(Y[Z_indice,])) )
S = s$d
- U = s$u
- V = s$v
#Set m-rank(r) singular values to zero, and recompose
#best rank(r) approximation of the initial product
if(rank[r] < length(S))
S[(rank[r]+1):length(S)] = 0
- phi[,,r] = U %*% diag(S) %*% t(V) %*% Rho[,,r]
+ phi[,,r] = s$u %*% diag(S) %*% t(s$v) %*% Rho[,,r]
}
-
+
#Etape E et calcul de LLF
sumLogLLF2 = 0
- for(i in 1:n){
+ for(i in seq_len(n))
+ {
sumLLF1 = 0
maxLogGamIR = -Inf
- for(r in 1:k){
+ for (r in seq_len(k))
+ {
dotProduct = tcrossprod(Y[i,]%*%Rho[,,r]-X[i,]%*%phi[,,r])
logGamIR = log(Pi[r]) + log(det(Rho[,,r])) - 0.5*dotProduct
#Z[i] = index of max (gam[i,])
- if(logGamIR > maxLogGamIR){
+ if(logGamIR > maxLogGamIR)
+ {
Z[i] = r
maxLogGamIR = logGamIR
}
- sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2)
+ sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2)
}
sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
}
LLF = -1/n * sumLogLLF2
-
+
#update distance parameter to check algorithm convergence (delta(phi, Phi))
- deltaPhi = c(deltaPhi, max(max(max((abs(phi-Phi))/(1+abs(phi))))) )
- if(length(deltaPhi) > deltaPhiBufferSize){
- l_1 = c(2:length(deltaPhi))
- deltaPhi = deltaPhi[l_1]
- }
+ deltaPhi = c( deltaPhi, max( (abs(phi-Phi)) / (1+abs(phi)) ) ) #TODO: explain?
+ if (length(deltaPhi) > deltaPhiBufferSize)
+ deltaPhi = deltaPhi[2:length(deltaPhi)]
sumDeltaPhi = sum(abs(deltaPhi))
-
+
#update other local variables
Phi = phi
ite = ite+1
}
- return(list(phi=phi, LLF=LLF))
+ return(list("phi"=phi, "LLF"=LLF))
}
+++ /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
-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
-checkOutput = function(varName, array, refArray, tol)
-{
- print(paste("Checking ",varName,sep=""))
- maxError = max(abs(array - refArray))
- if(maxError >= tol){
- print(paste("Inaccuracy: max(abs(error)) = ",maxError," >= ",tol,sep=""))
- } else{
- print("OK")
- }
-}
-constructionModelesLassoRank = function(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax){
+constructionModelesLassoRank = function(pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
+{
#get matrix sizes
n = dim(X)[1]
p = dim(X)[2]
m = dim(rho)[2]
k = dim(rho)[3]
L = dim(A1)[2]
-
+
+ # On cherche les rangs possiblement intéressants
deltaRank = rangmax - rangmin + 1
Size = deltaRank^k
- Rank = matrix(0, Size, k)
-# for(r in 1:k) {
-# Rank[,r] = rangmin + <--- #FIXME:
-# }
-
+ Rank = matrix(0, nrow=Size, ncol=k)
+ for(r in 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 :
+ # ça remplit la colonne
+ # Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre,
+ # et on fait ça (rangmax-rangmin)^2 fois
+ # ...
+ # Dans la dernière, on répète chaque chiffre une fois,
+ # et on fait ça (rangmin-rangmax)^(k-1) fois.
+ Rank[,r] = rangmin + rep(0:(deltaRank-1), deltaRank^(r-1), each=deltaRank^(k-r))
+ }
+
+ # output parameters
phi = array(0, dim=c(p,m,k,L*Size))
llh = matrix(0, L*Size, 2) #log-likelihood
- for(lambdaIndex in 1:L){
- #on ne garde que les colonnes actives
- #active sera l'ensemble des variables informatives
- active = A1[, lambdaIndex]
- active[active==0] = c()
- if(length(active)>0){
- for(j in 1:Size){
- EMG_rank = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi, X[, active], Y, tau, Rank[j,])
- phiLambda = EMG_rank$phi
- LLF = EMG_rank$LLF
- llh[(lambdaIndex-1)*Size+j,] = c(LLF, sum(Rank[j,]^(length(active)- Rank[j,]+m)))
- phi[active,,,(lambdaIndex-1)*Size+j] = phiLambda
+ for(lambdaIndex in 1:L)
+ {
+ # on ne garde que les colonnes actives
+ # 'active' sera l'ensemble des variables informatives
+ active = A1[,lambdaIndex]
+ active = active[-(active==0)]
+ if (length(active) > 0)
+ {
+ for (j in 1:Size)
+ {
+ res = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi,
+ X[,active], Y, tau, Rank[j,])
+ llh[(lambdaIndex-1)*Size+j,] =
+ c( res$LLF, sum(Rank[j,] * (length(active)- Rank[j,] + m)) )
+ phi[active,,,(lambdaIndex-1)*Size+j] = res$phi
}
}
}
- return(list(phi=phi, llh = llh))
+ return (list(phi=phi, llh = llh))
}
+++ /dev/null
-function[phi,llh] = 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);
- llh = 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,:));
- llh((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
-covariance = function(p,a)
-{
- A = matrix(a, p,p)
- for (i in 1:p){
- A[i,] = A[i,]^abs(i-(1:p))
- }
- return (A)
-}
+++ /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
#include "constructionModelesLassoRank.h"
#include "test_utils.h"
#include <stdlib.h>
+#include <math.h>
int main(int argc, char** argv)
{
////////////
// INPUTS //
- Real* Pi = readArray_real("Pi");
- Real* Rho = readArray_real("Rho");
+ Real* pi = readArray_real("pi");
+ Real* rho = readArray_real("rho");
int mini = read_int("mini");
int maxi = read_int("maxi");
Real* X = readArray_real("X");
/////////////////////////////////////////
// Call to constructionModelesLassoMLE //
constructionModelesLassoRank_core(
- Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
+ pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
phi,llh,
n,p,m,k,L);
/////////////////////////////////////////
- free(Rho);
- free(Pi);
+ free(rho);
+ free(pi);
free(X);
free(Y);
free(A1);
projects / valse.git / commitdiff