Commit | Line | Data |
---|---|---|
228ee602 | 1 | #' EMGLLF |
2 | #' | |
3 | #' Description de EMGLLF | |
4 | #' | |
5 | #' @param phiInit an initialization for phi | |
6 | #' @param rhoInit an initialization for rho | |
7 | #' @param piInit an initialization for pi | |
8 | #' @param gamInit initialization for the a posteriori probabilities | |
9 | #' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10 | |
10 | #' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100 | |
11 | #' @param gamma integer for the power in the penaly, by default = 1 | |
12 | #' @param lambda regularization parameter in the Lasso estimation | |
13 | #' @param X matrix of covariates (of size n*p) | |
14 | #' @param Y matrix of responses (of size n*m) | |
15 | #' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4 | |
16 | #' | |
17 | #' @return A list ... phi,rho,pi,LLF,S,affec: | |
18 | #' phi : parametre de moyenne renormalisé, calculé par l'EM | |
19 | #' rho : parametre de variance renormalisé, calculé par l'EM | |
20 | #' pi : parametre des proportions renormalisé, calculé par l'EM | |
21 | #' LLF : log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres | |
17b9fa5f BA |
22 | #' S : ... |
23 | #' affec : ... | |
228ee602 | 24 | #' |
25 | #' @export | |
26 | EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, | |
27 | X, Y, eps, fast) | |
28 | { | |
29 | if (!fast) | |
30 | { | |
31 | # Function in R | |
32 | return(.EMGLLF_R(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, | |
33 | X, Y, eps)) | |
34 | } | |
35 | ||
36 | # Function in C | |
37 | n <- nrow(X) #nombre d'echantillons | |
38 | p <- ncol(X) #nombre de covariables | |
39 | m <- ncol(Y) #taille de Y (multivarié) | |
40 | k <- length(piInit) #nombre de composantes dans le mélange | |
41 | .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, | |
42 | X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k), | |
43 | LLF = double(maxi), S = double(p * m * k), affec = integer(n), n, p, m, k, | |
44 | PACKAGE = "valse") | |
17b9fa5f | 45 | list(phi = phi, rho = rho, pi = pi, llh = llh, S = S, affec=affec) |
228ee602 | 46 | } |
47 | ||
48 | # R version - slow but easy to read | |
49 | .EMGLLF_R <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, | |
50 | X, Y, eps) | |
51 | { | |
52 | # Matrix dimensions | |
53 | n <- nrow(X) | |
54 | p <- ncol(X) | |
55 | m <- ncol(Y) | |
56 | k <- length(piInit) | |
57 | ||
58 | # Adjustments required when p==1 or m==1 (var.sel. or output dim 1) | |
59 | if (p==1 || m==1) | |
60 | phiInit <- array(phiInit, dim=c(p,m,k)) | |
61 | if (m==1) | |
62 | rhoInit <- array(rhoInit, dim=c(m,m,k)) | |
63 | ||
64 | # Outputs | |
65 | phi <- phiInit | |
66 | rho <- rhoInit | |
67 | pi <- piInit | |
68 | llh <- -Inf | |
69 | S <- array(0, dim = c(p, m, k)) | |
70 | ||
71 | # Algorithm variables | |
72 | gam <- gamInit | |
73 | Gram2 <- array(0, dim = c(p, p, k)) | |
74 | ps2 <- array(0, dim = c(p, m, k)) | |
75 | X2 <- array(0, dim = c(n, p, k)) | |
76 | Y2 <- array(0, dim = c(n, m, k)) | |
228ee602 | 77 | |
78 | for (ite in 1:maxi) | |
79 | { | |
80 | # Remember last pi,rho,phi values for exit condition in the end of loop | |
81 | Phi <- phi | |
82 | Rho <- rho | |
83 | Pi <- pi | |
84 | ||
85 | # Computations associated to X and Y | |
86 | for (r in 1:k) | |
87 | { | |
88 | for (mm in 1:m) | |
89 | Y2[, mm, r] <- sqrt(gam[, r]) * Y[, mm] | |
90 | for (i in 1:n) | |
91 | X2[i, , r] <- sqrt(gam[i, r]) * X[i, ] | |
92 | for (mm in 1:m) | |
93 | ps2[, mm, r] <- crossprod(X2[, , r], Y2[, mm, r]) | |
94 | for (j in 1:p) | |
95 | { | |
96 | for (s in 1:p) | |
97 | Gram2[j, s, r] <- crossprod(X2[, j, r], X2[, s, r]) | |
98 | } | |
99 | } | |
100 | ||
101 | ## M step | |
102 | ||
103 | # For pi | |
104 | b <- sapply(1:k, function(r) sum(abs(phi[, , r]))) | |
105 | gam2 <- colSums(gam) | |
106 | a <- sum(gam %*% log(pi)) | |
107 | ||
108 | # While the proportions are nonpositive | |
109 | kk <- 0 | |
110 | pi2AllPositive <- FALSE | |
111 | while (!pi2AllPositive) | |
112 | { | |
113 | pi2 <- pi + 0.1^kk * ((1/n) * gam2 - pi) | |
114 | pi2AllPositive <- all(pi2 >= 0) | |
115 | kk <- kk + 1 | |
116 | } | |
117 | ||
118 | # t(m) is the largest value in the grid O.1^k such that it is nonincreasing | |
119 | while (kk < 1000 && -a/n + lambda * sum(pi^gamma * b) < | |
120 | # na.rm=TRUE to handle 0*log(0) | |
121 | -sum(gam2 * log(pi2), na.rm=TRUE)/n + lambda * sum(pi2^gamma * b)) | |
122 | { | |
123 | pi2 <- pi + 0.1^kk * (1/n * gam2 - pi) | |
124 | kk <- kk + 1 | |
125 | } | |
126 | t <- 0.1^kk | |
127 | pi <- (pi + t * (pi2 - pi))/sum(pi + t * (pi2 - pi)) | |
128 | ||
129 | # For phi and rho | |
130 | for (r in 1:k) | |
131 | { | |
132 | for (mm in 1:m) | |
133 | { | |
134 | ps <- 0 | |
135 | for (i in 1:n) | |
136 | ps <- ps + Y2[i, mm, r] * sum(X2[i, , r] * phi[, mm, r]) | |
137 | nY2 <- sum(Y2[, mm, r]^2) | |
138 | rho[mm, mm, r] <- (ps + sqrt(ps^2 + 4 * nY2 * gam2[r]))/(2 * nY2) | |
139 | } | |
140 | } | |
141 | ||
142 | for (r in 1:k) | |
143 | { | |
144 | for (j in 1:p) | |
145 | { | |
146 | for (mm in 1:m) | |
147 | { | |
148 | S[j, mm, r] <- -rho[mm, mm, r] * ps2[j, mm, r] + | |
149 | sum(phi[-j, mm, r] * Gram2[j, -j, r]) | |
150 | if (abs(S[j, mm, r]) <= n * lambda * (pi[r]^gamma)) { | |
151 | phi[j, mm, r] <- 0 | |
152 | } else if (S[j, mm, r] > n * lambda * (pi[r]^gamma)) { | |
153 | phi[j, mm, r] <- (n * lambda * (pi[r]^gamma) - S[j, mm, r])/Gram2[j, j, r] | |
154 | } else { | |
155 | phi[j, mm, r] <- -(n * lambda * (pi[r]^gamma) + S[j, mm, r])/Gram2[j, j, r] | |
156 | } | |
157 | } | |
158 | } | |
159 | } | |
160 | ||
161 | ## E step | |
162 | ||
163 | # Precompute det(rho[,,r]) for r in 1...k | |
164 | detRho <- sapply(1:k, function(r) gdet(rho[, , r])) | |
165 | sumLogLLH <- 0 | |
166 | for (i in 1:n) | |
167 | { | |
168 | # Update gam[,]; use log to avoid numerical problems | |
169 | logGam <- sapply(1:k, function(r) { | |
170 | log(pi[r]) + log(detRho[r]) - 0.5 * | |
171 | sum((Y[i, ] %*% rho[, , r] - X[i, ] %*% phi[, , r])^2) | |
172 | }) | |
173 | ||
174 | logGam <- logGam - max(logGam) #adjust without changing proportions | |
175 | gam[i, ] <- exp(logGam) | |
176 | norm_fact <- sum(gam[i, ]) | |
177 | gam[i, ] <- gam[i, ] / norm_fact | |
178 | sumLogLLH <- sumLogLLH + log(norm_fact) - log((2 * base::pi)^(m/2)) | |
179 | } | |
180 | ||
181 | sumPen <- sum(pi^gamma * b) | |
182 | last_llh <- llh | |
183 | llh <- -sumLogLLH/n #+ lambda * sumPen | |
184 | dist <- ifelse(ite == 1, llh, (llh - last_llh)/(1 + abs(llh))) | |
185 | Dist1 <- max((abs(phi - Phi))/(1 + abs(phi))) | |
186 | Dist2 <- max((abs(rho - Rho))/(1 + abs(rho))) | |
187 | Dist3 <- max((abs(pi - Pi))/(1 + abs(Pi))) | |
188 | dist2 <- max(Dist1, Dist2, Dist3) | |
189 | ||
190 | if (ite >= mini && (dist >= eps || dist2 >= sqrt(eps))) | |
191 | break | |
192 | } | |
193 | ||
17b9fa5f BA |
194 | affec = apply(gam, 1, which.max) |
195 | list(phi = phi, rho = rho, pi = pi, llh = llh, S = S, affec=affec) | |
228ee602 | 196 | } |