Nothing
init_est_para.mfa <- function(Y, g, q, start, sigma_type = "common",
D_type = "common", ...) {
p <- ncol(Y)
n <- nrow(Y)
B <- array(NA, c(p, q, g))
pivec <- array(NA, c(g, 1))
mu <- array(NA, c(p, g))
if (D_type == 'common') {
D <- diag(diag(cov(Y)))
Di.sqrt <- diag(sqrt(diag(D)))
inv.Di.sqrt <- diag(1 / diag(Di.sqrt))
} else if (D_type == 'unique') {
D <- array(NA, c(p, p, g))
}
for(i in 1 : g) {
indices <- which(start == i)
pivec[i] <- length(indices) / n
mu[, i] <- apply(Y[indices, ], 2, mean)
Si <- cov(Y[indices, ])
if (D_type == 'unique'){
D[,, i] <- diag(diag(Si))
Di.sqrt <- diag(sqrt(diag(D[,, i])))
inv.Di.sqrt <- diag(1 / diag(Di.sqrt))
}
eig.list <- try(eigen(inv.Di.sqrt %*% Si %*% inv.Di.sqrt), TRUE)
if (is(eig.list, "try-error"))
break
H <- eig.list$vectors
sort.lambda <- sort(eig.list$values, decreasing = TRUE,
index.return = TRUE)
lambda <- sort.lambda$x
ix.lambda <- sort.lambda$ix
if (q == p) {
sigma2 <- 0
} else {
sigma2 <- mean(lambda[(q + 1) : p])
}
if (q == 1) {
B[,, i] <- Di.sqrt %*% H[, ix.lambda[1 : q]] %*%
diag((lambda[1 : q] - sigma2), q)
} else {
B[,, i] <- Di.sqrt %*% H[, ix.lambda[1 : q]] %*%
diag((lambda[1 : q] - sigma2))
}
}
if (sigma_type == 'common')
B <- apply(B, c(1, 2), sum) / g
model <- list(g = g, q = q, pivec = pivec, B = B, mu = mu, D = D,
sigma_type = sigma_type, D_type = D_type)
class(model) <- "mfa"
return(model)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.