228ee602 |
1 | #' initialization of the EM algorithm |
2 | #' |
3 | #' @param k number of components |
4 | #' @param X matrix of covariates (of size n*p) |
5 | #' @param Y matrix of responses (of size n*m) |
6 | #' |
7 | #' @return a list with phiInit, rhoInit, piInit, gamInit |
8 | #' @export |
9 | #' @importFrom methods new |
10 | #' @importFrom stats cutree dist hclust runif |
11 | initSmallEM <- function(k, X, Y, fast) |
12 | { |
13 | n <- nrow(X) |
14 | p <- ncol(X) |
15 | m <- ncol(Y) |
16 | nIte <- 20 |
17 | Zinit1 <- array(0, dim = c(n, nIte)) |
18 | betaInit1 <- array(0, dim = c(p, m, k, nIte)) |
19 | sigmaInit1 <- array(0, dim = c(m, m, k, nIte)) |
20 | phiInit1 <- array(0, dim = c(p, m, k, nIte)) |
21 | rhoInit1 <- array(0, dim = c(m, m, k, nIte)) |
22 | Gam <- matrix(0, n, k) |
23 | piInit1 <- matrix(0, nIte, k) |
24 | gamInit1 <- array(0, dim = c(n, k, nIte)) |
25 | LLFinit1 <- list() |
26 | |
27 | # require(MASS) #Moore-Penrose generalized inverse of matrix |
28 | for (repet in 1:nIte) |
29 | { |
30 | distance_clus <- dist(cbind(X, Y)) |
31 | tree_hier <- hclust(distance_clus) |
32 | Zinit1[, repet] <- cutree(tree_hier, k) |
33 | |
34 | for (r in 1:k) |
35 | { |
36 | Z <- Zinit1[, repet] |
37 | Z_indice <- seq_len(n)[Z == r] #renvoit les indices où Z==r |
38 | if (length(Z_indice) == 1) { |
39 | betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*% |
40 | crossprod(t(X[Z_indice, ]), Y[Z_indice, ]) |
41 | } else { |
42 | betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*% |
43 | crossprod(X[Z_indice, ], Y[Z_indice, ]) |
44 | } |
45 | sigmaInit1[, , r, repet] <- diag(m) |
46 | phiInit1[, , r, repet] <- betaInit1[, , r, repet] #/ sigmaInit1[,,r,repet] |
47 | rhoInit1[, , r, repet] <- solve(sigmaInit1[, , r, repet]) |
48 | piInit1[repet, r] <- mean(Z == r) |
49 | } |
50 | |
51 | for (i in 1:n) |
52 | { |
53 | for (r in 1:k) |
54 | { |
55 | dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet] |
56 | - X[i, ] %*% phiInit1[, , r, repet]) |
57 | Gam[i, r] <- piInit1[repet, r] * |
58 | gdet(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct) |
59 | } |
60 | sumGamI <- sum(Gam[i, ]) |
61 | gamInit1[i, , repet] <- Gam[i, ]/sumGamI |
62 | } |
63 | |
64 | miniInit <- 10 |
65 | maxiInit <- 11 |
66 | |
67 | init_EMG <- EMGLLF(phiInit1[, , , repet], rhoInit1[, , , repet], piInit1[repet, ], |
68 | gamInit1[, , repet], miniInit, maxiInit, gamma = 1, lambda = 0, X, Y, |
69 | eps = 1e-04, fast) |
70 | LLFinit1[[repet]] <- init_EMG$llh |
71 | } |
72 | b <- which.min(LLFinit1) |
73 | phiInit <- phiInit1[, , , b] |
74 | rhoInit <- rhoInit1[, , , b] |
75 | piInit <- piInit1[b, ] |
76 | gamInit <- gamInit1[, , b] |
77 | |
78 | return(list(phiInit = phiInit, rhoInit = rhoInit, piInit = piInit, gamInit = gamInit)) |
79 | } |