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