Nothing
fem.mstep <- function(Y,U,T,model,method){
# 12 different submodels: [DkBk] ... [AkjBk]
# Initialization
Y = as.matrix(Y)
n = nrow(Y)
p = ncol(Y)
K = ncol(T)
d = min(p-1,(K-1))
U = as.matrix(U)
mu = matrix(NA,K,d)
m = matrix(NA,K,p)
prop = rep(c(NA),1,K)
D = array(0,c(K,d,d))
b = rep(NA,K)
# Projection
X = as.matrix(Y) %*% as.matrix(U)
# Estimation
test = 0
bk = 0
for (k in 1:K){
nk = sum(T[,k])
if (nk ==0) stop("some classes become empty\n",call.=FALSE)
# Prior Probability
prop[k] = nk / n
# Mean in the latent space
mu[k,] = colSums((T[,k]*matrix(1,n,d))* X) / nk
# Observed space
m[k,] = colSums(T[,k]*matrix(1,n,p)* Y) / nk
YY = as.matrix(Y - t(m[k,]*matrix(1,p,n)))
if (nk < 1) denk = 1 else denk = (nk-1)
YYk = T[,k] * matrix(1,n,p) * YY
# Estimation of Delta k amongst 12 submodels
if (model %in% c('DkBk', 'DkB', 'DBk', 'DB')) {
D[k,,] = crossprod(YYk%*%U, YY%*%U) / denk
}
else if (model %in% c('AkjBk', 'AkjB', 'AjBk', 'AjB')) {
D[k,,] = diag(diag(crossprod(YYk %*%U, YY%*%U) / denk), d)
}
else if (model %in% c('AkBk', 'AkB', 'ABk', 'AB')) {
D[k,,] = diag(rep(Trace(crossprod(YYk %*%U, YY%*%U) / denk)/d,d), d)
}
# compute b[k]
if (model %in% c('DkB','DB','AkjB', 'AjB', 'AkB','AB')) {
bk = bk + (1 / n) * (sum(YYk*YY) - Trace(D[k,,]))
} else {
b[k] = (sum(YYk*YY) - Trace(D[k,,])) / (denk * (p-d))
}
}
# Compute Sigma after loop for homoscedastic models
Sigma = matrix(0, d, d)
if (model %in% c('DBk', 'DB', 'AjBk', 'AjB', 'AB', 'ABk')) {
for (k in 1:K) {
nk = sum(T[,k])
if (nk < 1) denk = 1 else denk = nk
Sigma = Sigma + (denk / n) * D[k,,]
}
for (k in 1:K) {
D[k,,] = Sigma
}
}
# compute beta for '*B' models
if (model %in% c('DkB','DB','AkjB','AkB','AjB','AB')) {
bk = bk / (p-d)
b = rep(bk, K)
}
# avoid numerical
b[b<=0] = 1e-3
for (k in 1:K) if (Trace(D[k,,]<1e-3)!=0) test = test+1
prms = list(K=K,p=p,mean=mu,my=m,prop=prop,D=D,b=b,model=model,method=method,V=U,test=test)
}
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.