1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
mvnormalmix.init = function (x, lambda = NULL, mu = NULL, sigma = NULL, k = 2, arbmean = TRUE, arbvar = TRUE)
{
n <- nrow(x)
p <- ncol(x)
y <- apply(x, 1, mean)
x <- x[order(y), ]
x.bin <- list()
for (j in 1:k) {
x.bin[[j]] <- x[max(1, floor((j - 1) * n/k)):ceiling(j *
n/k), ]
}
if (is.null(sigma)) {
if (arbvar) {
sigma.hyp = lapply(1:k, function(i) (apply(x.bin[[i]],
2, var))^-1)
sigma = lapply(1:k, function(i) diag(1/rexp(p, rate = sigma.hyp[[i]])))
}
else {
sigma.hyp = apply(sapply(1:k, function(i) (apply(x.bin[[i]],
2, var))^-1), 2, mean)
sigma = diag(1/rexp(p, rate = sigma.hyp))
}
}
if (is.null(sigma) == FALSE && arbvar == TRUE) {
k = length(sigma)
}
if (is.null(mu)) {
mu.hyp <- lapply(1:k, function(i) apply(x.bin[[i]], 2,
mean))
if (arbvar) {
mu <- lapply(1:k, function(i) as.vector(rmvnorm(1,
mu = as.vector(mu.hyp[[i]]), sigma = as.matrix(sigma[[i]]))))
}
else mu <- lapply(1:k, function(i) as.vector(rmvnorm(1,
mu = as.vector(mu.hyp[[i]]), sigma = as.matrix(sigma))))
if (arbmean==FALSE) {
mu <- apply(sapply(mu,as.vector),1,mean)
# mu <- lapply(1:k, function(i) mu)
}
}
if (is.null(mu) == FALSE && arbmean == TRUE){
k=length(mu)
}
if (is.null(lambda)) {
lambda <- runif(k)
lambda <- lambda/sum(lambda)
}
else k <- length(lambda)
list(lambda = lambda, mu = mu, sigma = sigma, k = k)
}
|
