3f58ea3544b52ba2a050455eee590c2370fa75f8
[valse.git] / R / selectVariables.R
1 #' selectVaribles
2 #' It is a function which construct, for a given lambda, the sets of
3 #' relevant variables and irrelevant variables.
4 #'
5 #' @param phiInit an initial estimator for phi (size: p*m*k)
6 #' @param rhoInit an initial estimator for rho (size: m*m*k)
7 #' @param piInit an initial estimator for pi (size : k)
8 #' @param gamInit an initial estimator for gamma
9 #' @param mini minimum number of iterations in EM algorithm
10 #' @param maxi maximum number of iterations in EM algorithm
11 #' @param gamma power in the penalty
12 #' @param glambda grid of regularization parameters
13 #' @param X matrix of regressors
14 #' @param Y matrix of responses
15 #' @param thres threshold to consider a coefficient to be equal to 0
16 #' @param tau threshold to say that EM algorithm has converged
17 #'
18 #' @return TODO
19 #'
20 #' @examples TODO
21 #'
22 #' @export
23 selectVariables <- function(phiInit,rhoInit,piInit,gamInit,
24 mini,maxi,gamma,glambda,X,Y,thres,tau)
25 {
26 dimphi <- dim(phiInit)
27 p <- dimPhi[1]
28 m <- dimPhi[2]
29 k <- dimPhi[3]
30 L <- length(glambda);
31 A1 <- array(0, dim <- c(p,m+1,L))
32 A2 <- array(0, dim <- c(p,m+1,L))
33 Rho <- array(0, dim <- c(m,m,k,L))
34 Pi <- array(0, dim <- c(k,L));
35
36 # For every lambda in gridLambda, comutation of the coefficients
37 for (lambdaIndex in c(1:L))
38 {
39 Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,
40 gamma,glambda[lambdaIndex],X,Y,tau);
41 phi <- Res$phi
42 rho <- Res$rho
43 pi <- Res$pi
44
45 # If a coefficient is larger than the threshold, we keep it
46 selectedVariables <- array(0, dim = c(p,m))
47 discardedVariables <- array(0, dim = c(p,m))
48 atLeastOneSelectedVariable <- false
49 for (j in c(1:p))
50 {
51 cpt <- 1
52 cpt2 <-1
53 for (mm in c(1:m))
54 {
55 if (max(abs(phi[j,mm,])) > thres)
56 {
57 selectedVariables[j,cpt] <- mm
58 cpt <- cpt+1
59 atLeastOneSelectedVariable <- true
60 } else
61 {
62 discardedVariables[j,cpt2] <- mm
63 cpt2 <- cpt2+1
64 }
65 }
66 }
67
68 # If no coefficients have been selected, we provide the zero matrix
69 # We delete zero coefficients: vec = indices of zero values
70 if (atLeastOneSelectedVariable)
71 {
72 vec <- c()
73 for (j in c(1:p))
74 {
75 if (selectedVariables(j,1) != 0)
76 vec <- c(vec,j)
77 # Else ( NOTE: [auder] else ?! TODO: explain? )
78 # we provide the indices of relevant coefficients
79 A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec)))
80 A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,]
81 A2[,1,lambdaIndex] <- 1:p
82 A2[,2:(m+1),lambdaIndex] <- discardedVariables
83 Rho[,,,lambdaIndex] <- rho
84 Pi[,lambdaIndex] <- pi
85 }
86 }
87 }
88
89 return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi))
90 }