Nothing
R/code.R
In pmcgd: pmcgd
Defines functions
class
.cov.wt.punzo
.m.step
.npar.model
.ncovpar
.model.type
.gpcm
.print.gpcm
.summary.gpcm
.plot.gpcm
.igpar
.igparv
.rgpar
.rwgpar
.igpark
.EM
.EMn
.EMv
.weights
.e.step
.loglik
.getall
.MAP
.combinewk
.msEEE
.sumSk.wt
.msEEV
.getA
.getOk
.getEkOk
.msVEV
.msVVV
.msEEI
.msVEI
.msEVI
.msVVI
.msEII
.msVII
.msVVE
.newD3.MM
.newD4.MM
.newD
.testval
.testgrad.D
.msEVE
.msVVE
.msVEE
.msEVV
.MCM
MS
Documented in
class
MS
#library("mixture")
#library("mnormt")
############################################
## Classification matrix given the groups ##
############################################
class <- function(groups,k){
n <- length(groups)
z <- array(0,c(n,k),dimnames=list(1:n,paste("comp",1:k,sep="")))
for(i in 1:n)
z[i,groups[i]] <- 1
return(z)
}
#################################################################
## Modified Weighted Covariance Matrix for contaminated models ##
#################################################################
.cov.wt.punzo <- function(x,wt,fact){
# fact contains the corrections factors due to the contamination
if (is.data.frame(x))
x <- as.matrix(x)
else if (!is.matrix(x))
stop("'x' must be a matrix or a data frame")
if (!all(is.finite(x)))
stop("'x' must contain finite values only")
p <- ncol(x)
n <- nrow(x)
mu <- array(colSums(wt*fact/sum(wt*fact)*x),c(p),dimnames=list(paste("X.",1:p,sep="")))
cov <- array(crossprod(sqrt(wt*fact/sum(wt))*(x-matrix(rep(mu,n),n,p,byrow=TRUE))),c(p,p),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep="")))
return(
list(
center = mu,
cov = cov
)
)
}
################################
## M-step of the EM algorithm ##
################################
.m.step <- function(data=NULL, covtype=NULL, w=NULL, fact=matrix(1,nrow(w),ncol(w)), v=1, D=NULL, mtol=NULL, mmax=NULL) {
G= ncol(w);
d= ncol(data);
Sk = array(0, c(d,d,G) )
gpar= list()
for (k in 1:G ) {
gpar[[k]] = list()
temp = .cov.wt.punzo(x=data, wt=w[,k], fact=fact[,k])
gpar[[k]]$mu = temp$center
if (!any(is.na( temp$cov))) gpar[[k]]$sigma = temp$cov
Sk[,,k] = temp$cov
}
gpar$pi = apply(w,2,mean)
temp = .model.type(modelname = covtype, Sk=Sk, ng=gpar$pi, D=D, mtol= mtol, mmax= mmax )
gpar$D = temp$D
for (k in 1:G ) {
gpar[[k]]$sigma = temp$sigma[,,k] #+ diag( 1-v , d, d)
gpar[[k]]$invSigma = temp$invSigma[,,k]
gpar[[k]]$logdet = temp$logdet[k]
}
return(gpar)
}
.npar.model <- function(modelname=NULL, p=NULL, G=NULL) {
val = numeric(3)
val[1] = G-1
val[2] = G*p
val[3] = .ncovpar(modelname= modelname, p=p, G=G)
val = sum(val)
return(val)
}
.ncovpar <- function(modelname=NULL, p=NULL, G=NULL) {
if (is.null(p)) stop("p is null")
if (is.null(G)) stop("G is null")
if (is.null(modelname)) stop("modelname is null")
if (modelname == "EII") npar = 1
else if (modelname == "VII") npar = G
else if (modelname == "EEI") npar = p
else if (modelname == "VEI") npar = p + G -1
else if (modelname == "EVI") npar = p*G - G +1
else if (modelname == "VVI") npar = p*G
else if (modelname == "EEE") npar = p*(p+1)/2
else if (modelname == "EEV") npar = G*p*(p+1)/2 - (G-1)*p
else if (modelname == "VEV") npar = G*p*(p+1)/2 - (G-1)*(p-1)
else if (modelname == "VVV") npar = G*p*(p+1)/2
else if (modelname == "EVE") npar = p*(p+1)/2 + (G-1)*(p-1)
else if (modelname == "VVE") npar = p*(p+1)/2 + (G-1)*p
else if (modelname == "VEE") npar = p*(p+1)/2 + (G-1)
else if (modelname == "EVV") npar = G*p*(p+1)/2 - (G-1)
else stop("modelname is not correctly defined")
return(npar)
}
.model.type <- function(modelname=NULL, Sk=NULL, ng=NULL, D=NULL, mtol=1e-10, mmax=10) {
if (is.null(modelname)) stop("modelname is null")
if (modelname == "EII") val = .msEII(Sk=Sk, ng=ng)
else if (modelname == "VII") val = .msVII(Sk=Sk, ng=ng)
else if (modelname == "EEI") val = .msEEI(Sk=Sk, ng=ng)
else if (modelname == "VEI") val = .msVEI(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "EVI") val = .msEVI(Sk=Sk, ng=ng)
else if (modelname == "VVI") val = .msVVI(Sk=Sk, ng=ng)
else if (modelname == "EEE") val = .msEEE(Sk=Sk, ng=ng)
else if (modelname == "EEV") val = .msEEV(Sk=Sk, ng=ng)
else if (modelname == "VEV") val = .msVEV(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "VVV") val = .msVVV(Sk=Sk, ng=ng)
else if (modelname == "EVE") val = .msEVE(Sk=Sk, ng=ng, D0=D, eplison= mtol, max.iter= mmax)
else if (modelname == "VVE") val = .msVVE(Sk=Sk, ng=ng, D0=D, eplison= mtol, max.iter= mmax)
else if (modelname == "VEE") val = .msVEE(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "EVV") val = .msEVV(Sk=Sk, ng=ng)
else stop("modelname or covtype is not correctly defined")
if (!is.list(val)) val = list(sigma=val)
return(val)
}
.gpcm <- function(data=NULL, G=1:3, mnames=NULL, start=0, label=NULL, veo=FALSE, nmax=1000, atol=1e-8, mtol=1e-8, mmax=10, pprogress=FALSE, pwarning=TRUE) {
if (is.null(data)) stop('Hey, we need some data, please! data is null')
if (!is.matrix(data)) stop('The data needs to be in matrix form')
if (!is.numeric(data)) stop('The data is required to be numeric')
if (nrow(data) == 1) stop('nrow(data) is equal to 1')
if (ncol(data) == 1) stop('ncol(data) is equal to 1; This function currently only works with multivariate data p > 1')
if (any(is.na(data))) stop('No NAs allowed.')
if (is.null(G)) stop('G is NULL')
G = as.integer(ceiling(G))
if (!is.integer(G)) stop('G is not a integer')
if ( any(G < 1)) stop('G is not a positive integer')
if (is.null(mnames) ) mnames = c("EII", "VII", "EEI", "VEI", "EVI", "VVI", "EEE", "EEV", "VEV", "VVV", "EVE", "VVE", "VEE", "EVV")
bic = array(0, dim= c(length(G), length(mnames), 3), dimnames=list(G, mnames, c('loglik', "npar", "BIC")) )
model = NULL; curBIC = Inf;
for (g in 1:length(G)) {
for (i in 1:length(mnames)) {
if ( pprogress ) print(c(G[g],mnames[i]))
if (veo | .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) < nrow(data)) {
a = try( { .EM(data=data, G=G[g], nmax=nmax, covtype=mnames[i], start=start, label=label, atol= atol, mtol=mtol, mmax=mmax ) }, pwarning)
if ( length(a) > 1) {
bic[g,i,1:2] = c(a$loglik[length(a$loglik)], .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = -2*bic[g,i,1] + bic[g,i,2]*log(nrow(data))
if ( bic[g,i,3] < curBIC) {
model = append(a,list(mtype=mnames[i]))
curBIC = bic[g,i,3]
}
} else {
bic[g,i,1:2] = c(NA, .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = NA
}
} else {
bic[g,i,1:2] = c(NA, .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = NA
}
}}
if ( is.null(start) ) startobject= "deterministic annealing with default values."
else if ( is.matrix(start)) startobject= "a user specified initialization matrix."
else if ( is.function(start)) startobject= "a user specified initialization function."
else if ( length(start) > 1 ) startobject= "deterministic annealing with user specified values"
else if ( start == 0) startobject= "k-means"
else if ( start > 0) startobject= paste( start, " random initializations", collapse="" )
bicModel = list(G=length(model$gpar$pi), covtype=model$mtype, bic=curBIC )
val = list( start=start, startobject= startobject, gpar=model$gpar, loglik=model$loglik, z=model$z, map=model$MAP, BIC=bic, bicModel= bicModel)
# val = list( model=model, BIC=bic, bicModel= bicModel)
class(val)<-"gpcm"
return(val)
}
.print.gpcm <-function(x, ...){
# a = object$BIC[,,1]
cat("The model choosen by applying the BIC criteria has ", x$bicModel$G, "component(s) and the ", x$bicModel$covtype, "covariance structure\n using ", x$startobject, "\n" )
#endPrint(a)
}
.summary.gpcm <- function(object, ...){
bicl = object$BIC[,,3]
cat("BIC for each model, number of components (rows), and covariance structure (columns).\n")
print.default( bicl)
}
.plot.gpcm <- function(x, ...) {
bicl = x$BIC[,,1]
g = dimnames(bicl)[[1]]
cov = dimnames(bicl)[[2]]
ncov = length(cov)
ng = length(g)
plot( g, bicl[,1], ylim=range(bicl,na.rm=TRUE), xlab="# components", ylab="BIC",type='l')
if (ncov >1) {
for (i in 2:ncov) lines( g, bicl[,i], col=i)
}
legend( "bottomright", legend=cov, col=1:ncov, lty=1)
}
.igpar = function(data=NULL, g=NULL, covtype=NULL, start=NULL, labels=NULL, mtol=NULL, mmax=NULL) {
if (is.null(start)) start = seq( 0.1, 1, length.out=5)
if (is.function(start) ) w = start(data=data, g=g, covtype=covtype)
else if ( !is.null(dim(start)) ) w = start
else if (length(start) > 1) w = .igparv(data=data, g=g, covtype=covtype, vseq0=start , mtol=mtol, mmax = mmax, labels=labels)
else if (start == 0) w = .igpark(data=data, g=g, covtype=covtype)
else if (start > 0) w = .rgpar(data=data, g=g, covtype=covtype, n=ceiling(start), labels=labels, mtol=mtol, mmax = mmax )
else stop(paste('Initialization method ', start, " does not compute!"))
if (!is.matrix(w)) stop("The zij initialization matrix is not a matrix")
if (any(w < 0)) stop("Some of the elements of the zij initialization matrix are less than zero")
# if ( any(apply(w,1,sum) !=1) ) stop("Some of the rows of the zij initialization matrix do not sume to 1")
if (nrow(w) < nrow(data) ) stop("The nrow(zij) of the initialization matrix is less than nrow(data)")
if (nrow(w) > nrow(data) ) stop("The nrow(zij) of the initialization matrix is greater than nrow(data)")
if (ncol(w) < g ) stop("The nrow(zij) of the initialization matrix is less than match g")
if (ncol(w) > g ) stop("The nrow(zij) of the initialization matrix is greater than match g")
w = .combinewk(w, label=labels)
gpar = .m.step(data=data, covtype=covtype, w=w, v=1, mtol=mtol, mmax = mmax)
return(gpar)
}
.igparv <- function(data=NULL, g=NULL,covtype=NULL, vseq0=NULL, labels=NULL, mtol=NULL, mmax = NULL) {
vseq0 = as.numeric(vseq0)
if (is.null(vseq0)) stop('The sequence for deterministic annealing is NULL')
if ( !all( vseq0 <=1 & vseq0 >=0) ) stop('The sequence for deterministic annealing must be between 0 and 1')
w = .rwgpar(data = data, g=g, covtype = covtype)
w = .combinewk(w, label=labels)
gpar = .m.step(data=data, covtype=covtype, w=w, v=1, mtol= mtol, mmax = mmax)
gpar = .EMv(data=data, gpar0 = gpar, G=g, vseq=vseq0, m=1, label=labels, covtype = covtype, mtol= mtol, mmax = mmax )$gpar
w = .weights(data=data, gpar= gpar)
return(w)
}
.rgpar <- function(data=NULL, g=NULL, covtype=NULL, n=1, labels=NULL, mtol=NULL, mmax = NULL) {
w.old = .rwgpar(data= data, g=g,covtype= covtype, labels=labels)
gpar.old = .m.step(data=data, covtype=covtype, w=w.old, v=1, mtol=mtol, mmax = mmax)
loglik.old = .loglik(data=data, gpar=gpar.old)
for (i in 1:n) {
w.new = .rwgpar(data= data, g=g,covtype= covtype, labels=labels)
gpar.new = .m.step(data=data, covtype=covtype, w=w.new, v=1, mtol=mtol, mmax = mmax)
loglik.new = .loglik(data, gpar.new)
if (loglik.new > loglik.old) {
w.old = w.old
gpar.old = gpar.new
loglik.old = loglik.new
}
}
return(w.old)
}
.rwgpar <- function(data=NULL, g=NULL,covtype=NULL, labels=NULL) {
w = matrix(rexp(nrow(data)*g),nrow=nrow(data),ncol=g)
w = matrix(t(apply(w,1, function(z) { z/sum(z)})),nrow=nrow(data),ncol=g)
w = .combinewk(w, label=labels)
return(w)
}
.igpark <- function(data=NULL, g=NULL,covtype=NULL) {
lw = kmeans(data, centers=g, iter.max=10)$cluster
w = .combinewk(matrix(0,nrow=nrow(data),ncol=g), label=lw)
return(w)
}
.EM <- function(data=NULL, gpar0=NULL, G=2, start=1, label=NULL, covtype=NULL, nmax=1000, atol=1e-8, mtol=1e-8, mmax=10 ) {
val = list()
if (is.null(gpar0)) val$gpar = .igpar(data=data, g=G, covtype=covtype, start=start, labels=label, mtol=mtol, mmax=mmax)
else val$gpar = gpar0
val$loglik = numeric(nmax)
val$loglik[1] = .loglik(data=data, gpar=val$gpar)
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[2] = .loglik(data=data, gpar=val$gpar)
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[3] = .loglik(data=data, gpar=val$gpar)
i = 3
while ( ( .getall(val$loglik[(i-2):i]) > atol) & (i < (nmax) ) ) {
i = i+1
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, D=val$gpar$D, mtol=mtol, mmax=mmax)
val$loglik[i] = .loglik(data=data, gpar=val$gpar)
}
val$loglik = val$loglik[1:i]
val$z = .e.step(data=data, gpar=val$gpar, labels=label)
val$MAP = .MAP(data=data, gpar=val$gpar, label=label)
return(val)
}
.EMn <- function(data=NULL, gpar0 = NULL, G=2, n=10, label =NULL, covtype="VVV", mtol=1e-8, mmax=10 ) {
val = list()
if (is.null(gpar0)) val$gpar = .igpar(data=data, g=G, covtype=covtype, mtol=mtol, mmax=mmax)
else val$gpar = gpar0
val$loglik = numeric(n)
for (i in 1:n) {
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[i] = .loglik(data=data, gpar=val$gpar)
}
# val$label = label
return(val)
}
.EMv <- function(data=NULL, gpar0=NULL, G=3, vseq=c(1,1), m=2, label=NULL, covtype="VVV", mtol=NULL, mmax = NULL ) {
val = list()
if (is.null(gpar0)) val$gpar = .rgpar(data=data, g=G, covtype=covtype, n=1 )
else val$gpar = gpar0
val$loglik = numeric(length(vseq)*m)
count = 1
for (i in 1:length(vseq)) { for (j in 1:m) {
tempw = .e.step(data=data, gpar=val$gpar, labels=label, v=vseq[i])
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, D=val$gpar$D, mtol= mtol, mmax = mmax)
val$loglik[i] = .loglik(data, val$gpar)
count = count + 1
#cat("iteration = ", count-1, "\t v= ", vseq[i], "\t loglik = ", val$loglik[count-1], "\t pi=", val$gpar$pi, "\n")
}}
return(val)
}
.weights <- function(data=NULL, gpar=NULL, v=1) {
d = ncol(data)
G = length(gpar$pi)
if (G > 1) {
zlog = matrix(0, nrow=nrow(data), ncol=length(gpar$pi))
for (k in 1:G ) zlog[,k] = -1/2*mahalanobis(x=data, center=gpar[[k]]$mu, cov=gpar[[k]]$invSigma, inverted=TRUE) -1/2*gpar[[k]]$logdet - d/2 *log(2*pi)
w = t(apply( zlog, 1, function(z,wt,v) {
x= exp( v*(z + log(wt)) )
x=x/sum(x);
return(x) }, wt=gpar$pi,v=v ))
# w = t(apply(zlog, 1, function(z,wt,v) {
# x=exp(v*(z + log(wt)) );
# if (sum(x) == 0) x= rep(1,length(x))
# x = x/sum(x)
# return( x )
# }, wt=gpar$pi,v=v ))
} else w = matrix(1,nrow=nrow(data), ncol=G)
return(w)
}
.e.step <- function(data=NULL, gpar=NULL, labels=NULL, v=1) {
w = .weights(data=data, gpar=gpar,v=v)
if (!is.null(labels)) w = .combinewk(weights=w, label= labels)
return(w)
}
.loglik = function(data, gpar) {
# output is a G x nrow(data) matrix
d = ncol(data)
G = length(gpar$pi)
zlog = matrix(0, nrow=nrow(data), ncol=G)
for (k in 1:G) zlog[,k] = -1/2*mahalanobis(x=data, center=gpar[[k]]$mu, cov=gpar[[k]]$invSigma, inverted=TRUE) -1/2*gpar[[k]]$logdet - d/2*(log(2)+log(pi))
w = apply( exp(zlog),1,function(z,wt) { sum(z*wt) } , wt=gpar$pi)
val = sum(log(w))
if( is.nan(val) ) val =NA
return(val)
}
.getall <- function(loglik) {
lm1 = loglik[3]
lm = loglik[2]
lm_1 = loglik[1]
am = (lm1 - lm)/(lm - lm_1)
lm1.Inf = lm + (lm1 - lm)/(1-am)
val = lm1.Inf - lm
if (is.nan(val) ) val =0
return( abs(val) )
}
.MAP <- function(data, gpar, label=NULL) {
w = .weights(data=data, gpar=gpar, v=1)
if (!is.null(label)) w = .combinewk(weights=w, label= label)
z = apply(w, 1, function(z) { z=(1:length(z))[z==max(z)]; return(z[1]) })
z = as.numeric(z)
return( z)
}
.combinewk <- function(weights=NULL, label=NULL) {
# known is a numeric with
# 0 if unknown group membership
# 1,2,3,.. for label of known group
if (!is.null(label)) { #stop('label is null')
label = as.integer(label)
if (any(!is.integer(label))) stop("Labels are not integers")
if (any(label < 0 )) stop("Labels can only be positive integers")
if (ncol(weights) < max(label) ) stop("Number of groups is less then the number groups given by labels")
if ( sum(label!=0) == nrow(weights) ) {
if (ncol(weights) > max(label) ) stop("Every observations has a label; Cannot fit more groups to the data then given the by the labels.")
}
kw = label !=0
for (j in 1:ncol(weights)) weights[kw,j] = (label == j)[kw]
}
return(weights)
}
.msEEE <- function(Sk=NULL, ng=NULL) {
# Sk is an array of with dim (p x p x G)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = W
logdetW = log(det(W))
invW = solve(W)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = W
val$invSigma[,,g] = invW
val$logdet[g] = logdetW
}
return(val)
}
.sumSk.wt <- function(Sk=NULL, wt=NULL, d=NULL, G=NULL) {
# Sum Sk over the groups with weights wt.
W = matrix(0, nrow=d, ncol=d )
for (g in 1:G) W = W + Sk[,,g]* wt[g]
return(W)
}
.msEEV <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
EWk = Sk
A = matrix(0,d,d)
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk[,,g] = eigen(Wk)$vectors
A = A + t(EWk[,,g]) %*% Wk %*% EWk[,,g]
}
lam = prod(diag(A))^(1/d)
A = A/lam
lam = lam/sum(ng)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam *( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val$invSigma[,,g] = 1/lam *( EWk[,,g] %*% diag(1/diag(A),d) %*% t(EWk[,,g]) )
val$logdet[g] = d*log(lam)
}
return(val)
}
.getA <- function(Ok=NULL, lam=NULL, G=NULL, d=NULL) {
A = matrix(0,d,d)
for (g in 1:G) A = A + Ok[,,g]/lam[g]
A = diag(A)
A = diag( A/prod(A)^(1/d) )
return( A )
}
.getOk <- function(Sk=NULL, ng=NULL, G=NULL) {
Ok = Sk
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk = eigen(Wk)$vectors
Ok[,,g] = t(EWk) %*% Wk %*% EWk
Ok[,,g] = diag(diag(Ok[,,g]))
}
return(Ok)
}
.getEkOk <- function(Sk=NULL, ng=NULL, G=NULL) {
Ok = Sk
EWk = Sk
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk[,,g] = eigen(Wk)$vectors
Ok[,,g] = t(EWk[,,g]) %*% Wk %*% EWk[,,g]
}
return(list(Ok=Ok,EWk=EWk))
}
.msVEV <- function(Sk=NULL, ng=NULL, eplison=1e-14, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
temp = .getEkOk(Sk=Sk, ng=ng, G=G)
Ok = temp$Ok
EWk = temp$EWk
lam = apply(Ok,3,function(z) { sum(diag(z)) } )/(ng*d)
A = .getA(Ok=Ok, lam=lam, d=d, G=G)
lam = apply(Ok,3,function(z, invA) { sum(diag(z*invA)) }, invA=diag(1/diag(A)) )/(ng*d)
conv = c( d*sum(ng*(1+log(lam))), Inf )
count = 1
while ( diff(conv)/conv[1] > eplison & count < max.iter) {
A = .getA(Ok=Ok, lam=lam, d=d, G=G)
lam = apply(Ok,3,function(z, invA) { sum(diag(z*invA)) }, invA=diag(1/diag(A)) )/(ng*d)
conv = c(d*sum(ng*(1+log(lam))), conv[1] )
count = count +1
}
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = lam[g] * ( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g] * ( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val$invSigma[,,g] = 1/lam[g] * ( EWk[,,g] %*% diag(1/diag(A),d) %*% t(EWk[,,g]) )
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msVVV <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = Sk[,,g]
val$invSigma[,,g] = solve(Sk[,,g])
val$logdet[g] = log(det(Sk[,,g]) )
}
return(val)
}
.msEEI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
B = diag(diag(W))
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = B
val$invSigma[,,g] = diag( 1/diag(B),d )
val$logdet[g] = sum(log( diag(B) ))
}
return(val)
}
.msVEI <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
lam = apply(Sk,3,function(z) { sum(diag(z)) } )/d
# W = sumSk.wt(Sk=Sk, wt=ng*lam, d=d, G=G)
W = .sumSk.wt(Sk=Sk, wt=ng/lam, d=d, G=G)
W = diag(W)
B = diag(W/prod(W)^(1/d))
lam = apply(Sk,3,function(z, invB) { sum(diag(z*invB)) }, invB=diag(1/diag(B) ) )/d
conv = c( d*sum(ng*(1+log(lam))), Inf )
count = 1
while ( abs(diff(conv)) > eplison & count < max.iter) {
# while ( count < max.iter) {
# W = sumSk.wt(Sk=Sk, wt=ng*lam, d=d, G=G)
W = .sumSk.wt(Sk=Sk, wt=ng/lam, d=d, G=G)
W = diag(W)
B = diag(W/prod(W)^(1/d))
lam = apply(Sk,3,function(z, invB) { sum(diag(z*invB)) }, invB=diag(1/diag(B) ) )/d
conv = c(d*sum(ng*(1+log(lam))), conv[1] )
count = count +1
}
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g]*B
val$invSigma[,,g] = diag( 1/diag(B) * 1/lam[g],d )
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msEVI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Bk = matrix(0,nrow=d, ncol=G)
for (g in 1:G) Bk[,g] = diag(Sk[,,g]*ng[g])
lam = apply(Bk, 2, prod)^(1/d)
Bk = sweep(Bk, 2, 1/lam, FUN="*")
lam = sum(lam)/sum(ng)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam * diag(Bk[,g],d)
val$invSigma[,,g] = 1/lam * diag(1/Bk[,g],d)
val$logdet[g] = d*log(lam)
}
return(val)
}
.msVVI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
# Bk = matrix(0,nrow=d, ncol=G)
# for (g in 1:G) Bk[,g] = diag(Sk[,,g]*ng[g])
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = diag(diag(Sk[,,g]),d)
val$invSigma[,,g] = diag(1/diag(Sk[,,g]),d)
val$logdet[g] = sum( log(diag(Sk[,,g])) )
}
return(val)
}
.msEII <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
lam = sum(diag(W))/(sum(ng)*d)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = diag(rep(lam,d), d)
val$invSigma[,,g] = diag(rep(1/lam,d),d)
val$logdet[g] = d*log(lam)
}
return(val)
}
.msVII <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
sumdiagSkg = sum(diag(Sk[,,g]))
val$sigma[,,g] = diag(rep(sumdiagSkg/d,d))
val$invSigma[,,g] = diag(rep( d/sumdiagSkg, d))
val$logdet[g] = d* log( sumdiagSkg ) - d*log(d)
}
return(val)
}
.msVVE <- function(Sk=NULL, ng=NULL) {
# Sk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
# lam = numeric(G)
# for (g in 1:G) lam[g] = sum(diag(Sk[,,g]))/d
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = diag(rep(sum(diag(Sk[,,g]))/d,d))
return(val)
}
#####################
### new.cov models
#####################
.newD3.MM <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
z = matrix(0,d,d)
lambda =0
for (g in 1:G) {
lambdak = max(eigen(Wk[,,g])$values)
z = z + diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g] -lambdak *(diag(1/Ak[,g])%*% t(D) )
}
z1 = svd(z)
Xk1 = (z1$v) %*% t(z1$u)
return( Xk1 )
}
.newD4.MM <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
z = matrix(0,d,d)
lambda =0
for (g in 1:G) {
lambdak = max(1/Ak[,g])
#z = z + diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g] -lambdak *(diag(1/Ak[,g])%*% t(D) )
z = z + Wk[,,g] %*% (D) %*% diag(1/Ak[,g]) -lambdak *(Wk[,,g]%*% (D) )
}
z1 = svd(z)
# Xk1 = (z1$v) %*% t(z1$u)
# return( t(Xk1) )
# OR
Xk1 = (z1$v) %*% t(z1$u)
return( t(Xk1) )
}
.newD <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
D6 = D
D6 = .newD3.MM(D=D6, d=d, G=G, Wk=Wk, Ak=Ak, tmax = 100)
D6 = .newD4.MM(D=D6, d=d, G=G, Wk=Wk, Ak=Ak, tmax = 100)
return(D6)
}
.testval <- function(Wk=NULL, Ak=NULL, D=NULL, G=NULL) {
z = numeric(G)
# for (g in 1:G) z[g] = sum(diag( D %*% diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g]))
for (g in 1:G) z[g] = sum(diag( t(D) %*% Wk[,,g] %*% D %*% diag(1/Ak[,g]) ))
return(sum(z))
}
.testgrad.D <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL) {
z = matrix(0,d,d)
for (g in 1:G) z = z + Wk[,,g] %*% D %*% diag(1/Ak[,g])
return(2*z)
}
.msEVE <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 10, D0=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1]; G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
# D = diag(rep(1,d))
D = D0
if (is.null(D)) D = diag(rep(1,d))
# D = t(eigen(W)$vectors)
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% (D) ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(0, 1, testval(Wk=Wk,Ak=Ak,D=D,G=G)))
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
#print(c(0, 2, testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c( .testval(Wk=Wk,Ak=Ak,D=D,G=G), Inf )
count = 1
while ( diff(conv)/abs(conv[1]) > eplison & count < max.iter) {
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% (D) ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(count, 0, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c(.testval(Wk=Wk,Ak=Ak,D=D,G=G), conv[1] )
count = count +1
}
lam = 0
for (g in 1:G) lam = lam + sum(diag( D %*% diag(1/Ak[,g])%*% t(D) %*% Wk[,,g] ))
lam = lam/(sum(ng)*d)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G), D=D )
for (g in 1:G) {
val$sigma[,,g] = D %*% diag(lam*Ak[,g])%*% t(D)
val$invSigma[,,g] = D %*% diag(1/lam*1/Ak[,g])%*% t(D)
val$logdet[g] = d*log( lam )
}
return(val)
}
.msVVE <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 10, D0=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1]; G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
D = D0
if (is.null(D)) D = diag(rep(1,d))
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% D ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(0, 1, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
#print(c(0, 2, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c( .testval(Wk=Wk,Ak=Ak,D=D,G=G), Inf )
count = 1
while ( diff(conv)/abs(conv[1]) > eplison & count < max.iter) {
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% D ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
conv = c(.testval(Wk=Wk,Ak=Ak,D=D,G=G), conv[1] )
count = count +1
}
#print(count)
lam = numeric(G)
for (g in 1:G) lam[g] =sum(diag( D %*% diag(1/Ak[,g])%*% t(D) %*% Sk[,,g] ))/d
#print(apply(Ak,2,prod))
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = (D %*% diag(lam[g]*Ak[,g])%*% t(D))
val$invSigma[,,g] = (D %*% diag(1/lam[g]*1/Ak[,g])%*% t(D))
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msVEE <- function(Sk=NULL, ng=NULL, eplison=1e-14, max.iter=100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
C = W/det(W)^(1/d)
invC = solve(C)
lam = apply(Sk,3,function(z, invC) { sum(diag(z %*% invC)) }, invC=invC)/d
val1 = sum(apply(Wk,3,function(z, invC) { sum(diag(z*invC)) }, invC= invC)/lam) +d*sum(ng*lam)
conv = c(val1, Inf )
count = 1
while ( diff(conv)/conv[1] > eplison & count < max.iter) {
#for (i in 1:max.iter) {
C = .sumSk.wt(Sk=Wk, wt=1/lam, d=d,G=G)
C = C/det(C)^(1/d)
invC = solve(C)
lam = apply(Sk,3,function(z, invC) { sum(diag(z %*% invC)) }, invC=invC)/d
val1 = sum(apply(Wk,3,function(z, invC) { sum(diag(z*invC)) }, invC= invC )/lam) +d*sum(ng*lam)
conv = c(val1, conv[1] )
count = count +1
}
#print(c(count,det(C), lam) )
invC = solve(C)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g]*C
val$invSigma[,,g] = 1/lam[g]*invC
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msEVV <- function(Sk=NULL, ng=NULL, eplison=1e-12, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Wk = Sk
Ck = Sk
lam = numeric(G)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
lam[g] = det(Wk[,,g])^(1/d)
}
Ck = sweep(Wk, 3, 1/lam, FUN="*")
lam = sum(lam)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam * Ck[,,g]
val$invSigma[,,g] = 1/lam * solve(Ck[,,g])
val$logdet[g] = d* log(lam )
}
return(val)
}
####################################
## Mixture of Contaminated Models ##
####################################
.MCM <- function(
X, # matrix of data
k, # number of groups
initialization="mclust", # initialization procedure: "random.soft", "random.hard", "manual", or "mclust" (which is nested)
modelname="VVV", # one of the 14 models of Celeaux & Govaert (1995)
alphacon=TRUE, # if TRUE, alpha is constrained to be >0.5
alphamin=rep(0.5,k), # minimum proportion of good data in each group
alphafix=FALSE, # if TRUE, the inner weights are not estimated
alpha=NULL, # vector of dimension k with proportion of good observations in each group
etacon=TRUE, # if TRUE, eta is constrained to be >1
etafix=FALSE, # if TRUE, eta is not estimated
eta=NULL, # inflation parameters
etamax=200, # maximum value of eta
start.z=NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual"
start.v=NULL, # (n x 2 x k)-array of soft or hard classification in each group: it is used only if initialization="manual"
start=0, # initialization for the package mixture
ind.label=NULL, # indexes of the labelled observations
label=NULL, # groups of the labelled observations
iter.max=1000, # maximum number of iterations in the EM-algorithm
threshold=1.0e-04 # stopping rule in the Aitken rule
#loglikplot=TRUE, # if TRUE, the log-likelihood values against the iterations are plotted
#plot=TRUE # if TRUE, the plot of the classified data is given
)
{
#call = match.call()
if (is.data.frame(X))
X <- as.matrix(X)
n <- nrow(X) # sample size
p <- ncol(X) # number of variables
if (is.null(X)) stop('Hey, we need some data, please! X is null')
if (!is.matrix(X)) stop('X needs to be in matrix form')
if (!is.numeric(X)) stop('X is required to be numeric')
if (n == 1) stop('nrow(X) is equal to 1')
if (p == 1) stop('ncol(X) is equal to 1; This function currently only works with multivariate data p > 1')
if (any(is.na(X))) stop('No NAs allowed.')
if (is.null(k)) stop('k is NULL')
k <- as.integer(ceiling(k))
if (!is.integer(k)) stop('k is not a integer')
if (any(k < 1)) stop('k is not a positive integer')
# for model-based classification
lab <- NULL
if(is.vector(label)){
nlab <- length(label)
nunlab <- n-nlab
lab <- numeric(n)
lab[ind.label] <- label
}
prior <- numeric(k) # proportion of each group
priorgood <- numeric(k) # proportion of good observations in each group
v <- array(0.5,c(n,2,k),dimnames=list(1:n,c("good","bad"),paste("group ",1:k,sep="")))
mu <- array(0,c(p,k),dimnames=list(paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
Sigma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
invSigma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
lambda <- array(0,c(k),dimnames = list(paste("group ",1:k,sep="")))
Delta <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
Gamma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
if(!etafix)
eta <- rep(2,k) # inflation parameters
correction <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep=""))) # factor which differentiates this model by mclust
PXgood <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep="")))
PXbad <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep="")))
# ------------------------- #
# posteriors initialization #
# ------------------------- #
if(initialization=="random.soft"){
z <- array(runif(n*k),c(n,k)) # soft posterior probabilities (no-normalized) (n x k)
z <- z/rowSums(z) # soft posterior probabilities (n x k)
for(j in 1:k){
v[,2,j] <- runif(n,0,0.2)
v[,1,j] <- 1-v[,2,j]
}
}
if(initialization=="random.hard"){
z <- t(rmultinom(n, size = 1, prob=rep(1/k,k))) # hard posterior probabilities (n x k)
for(j in 1:k)
v[,,j] <- t(rmultinom(n, size = 1, prob=c(0.9,0.1)))
}
if(initialization=="manual"){
z <- start.z # soft or hard classification
v <- start.v # soft or hard classification of good and bad observations
}
if(initialization=="mclust"){
mclustfit <- .gpcm(data=X, G=k, mnames=modelname, start=start, label=lab, pprogress=FALSE)
for(j in 1:k){
mu[,j] <- mclustfit$gpar[[j]]$mu
Sigma[,,j] <- mclustfit$gpar[[j]]$sigma
invSigma[,,j] <- mclustfit$gpar[[j]]$invSigma
#temp <- eigen(Sigma[,,j])
#Gamma[,,j] <- temp$vectors
#Delta[,,j] <- diag(temp$values)/prod(temp$values)^(1/p)
#lambda[j] <- prod(temp$values)^(1/p)
}
z <- mclustfit$z
group <- apply(z,1,which.max)
for(j in 1:k){
Probj <- (2*pi)^(-p/2)*det(Sigma[,,j])^(-1/2)*exp(-1/2*mahalanobis(x=X[which(group==j),], center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
#Probj <- dmnorm(x=X[which(group==j),], mean = mu[,j], varcov=Sigma[,,j])
wgood <- Probj/max(Probj)
v[which(group==j),1,j] <- (wgood)^0.2
v[which(group==j),2,j] <- 1-(wgood)^0.2
}
}
if(is.vector(label))
z[ind.label,] <- class(label,k=k)
# ------------ #
# EM algorithm #
# ------------ #
# Preliminary definition of convergence criterions
check <- 0
iteration <- 1
loglik <- NULL
aloglik <- NULL
aloglik <- c(0,0)
a <- NULL
a <- c(0,0)
while(check<1){
# ++++++ #
# M-step #
# ++++++ #
# ------- #
# Weights #
# ------- #
prior <- colMeans(z)
zv <- z*v[,1,]
zvcompl <- z*v[,2,]
if(!alphafix){
if(alphacon){
g <- function(alpha,j,z,v)
sum(z[,j]*(v[,1,j]*log(alpha)+v[,2,j]*log(1-alpha)))
for(j in 1:k)
priorgood[j] <- optimize(g,c(alphamin[j],1),maximum = TRUE,j=j,z=z,v=v)$maximum
}
if(!alphacon)
priorgood <- colSums(zv)/colSums(z)
}
if(alphafix)
priorgood <- alpha
priorbad <- 1-priorgood
# ---------- #
# mu & Sigma #
# ---------- #
correction <- v[,1,]+v[,2,]*matrix(rep(1/eta),n,k,byrow=TRUE)
fitM <- .m.step(data=X, covtype=modelname, w=z, fact=correction, v=1, mtol=1e-10, mmax=10)
for(j in 1:k){
mu[,j] <- fitM[[j]]$mu
Sigma[,,j] <- fitM[[j]]$sigma
invSigma[,,j] <- fitM[[j]]$invSigma
temp <- eigen(Sigma[,,j])
Gamma[,,j] <- temp$vectors
Delta[,,j] <- diag(temp$values)/prod(temp$values)^(1/p)
lambda[j] <- prod(temp$values)^(1/p)
}
# ------------------- #
# Inflation parameter #
# ------------------- #
if(!etafix){
if(etacon){
f <- function(eta,j,p,zvcompl,X,mu,invSigma)
sum(-(p/2)*zvcompl[,j]*log(eta)-(1/2)*zvcompl[,j]*(1/eta)*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
for(j in 1:k)
eta[j] <- optimize(f,c(1,etamax),maximum = TRUE,j=j,p=p,zvcompl=zvcompl,X=X,mu=mu,invSigma=invSigma)$maximum
}
if(!etacon){
for(j in 1:k)
eta[j] <- sum(zvcompl[,j]*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))/(p*sum(zvcompl[,j]))
}
}
# ------- #
# density #
# ------- #
zerostar <- 1e-323 # to avoid zero probabilities
for(j in 1:k){
PXgood[,j] <- (2*pi)^(-p/2)*(det(Sigma[,,j]))^(-1/2)*exp(-1/2*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
PXbad[,j] <- (2*pi)^(-p/2)*(det(eta[j]*Sigma[,,j]))^(-1/2)*exp(-1/2*mahalanobis(x=X, center=mu[,j], cov=1/eta[j]*invSigma[,,j], inverted=TRUE))
PXgood[,j] <- (PXgood[,j]<zerostar)*zerostar+(PXgood[,j]>=zerostar)*PXgood[,j]
PXbad[,j] <- (PXbad[,j]<zerostar)*zerostar+(PXbad[,j]>=zerostar)*PXbad[,j]
}
# ------------------------------------- #
# Global - Observed-data log-likelihood #
# ------------------------------------- #
# model-based clustering
if(!is.vector(label))
llvalues <- sum(log( rowSums(matrix(rep(prior,n),n,k,byrow=TRUE)*(matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad))))
# model-based classification
if(is.vector(label)){
llvalueslab <- z[ind.label,]*(log(matrix(rep(prior,nlab),nlab,k,byrow=T))+log( matrix(rep(priorgood,nlab),nlab,k,byrow=TRUE)*PXgood[ind.label,]+matrix(rep(priorbad,nlab),nlab,k,byrow=TRUE)*PXbad[ind.label,] ))
llvaluesunlab <- log(rowSums(matrix(rep(prior,nunlab),nunlab,k,byrow=TRUE)*(matrix(rep(priorgood,nunlab),nunlab,k,byrow=TRUE)*PXgood[-ind.label,]+matrix(rep(priorbad,nunlab),nunlab,k,byrow=TRUE)*PXbad[-ind.label,])))
llvalues <- sum(llvalueslab)+sum(llvaluesunlab)
}
loglik[iteration] <- llvalues
# ----------------------------------------------- #
# Aitkane's Acceleration-Based Stopping Criterion #
# ----------------------------------------------- #
if(iteration>2 & k > 1){
if(abs(loglik[iteration-1]-loglik[iteration-2])>0){
a[iteration-1] <- (loglik[iteration]-loglik[iteration-1])/(loglik[iteration-1]-loglik[iteration-2])
aloglik[iteration] <- loglik[iteration-1]+(1/(1-a[iteration-1])*(loglik[iteration]-loglik[iteration-1]))
if(abs(aloglik[iteration]-loglik[iteration])<threshold)
check <- 1
}
#if(abs(loglik[iteration-1]-loglik[iteration-2])==0)
else
check <- 1
}
if(iteration==iter.max | k==1) check <- 1
cat("*")
#cat("Iteration",iteration,"\n")
iteration <- iteration + 1
# ++++++ #
# E-Step #
# ++++++ #
z.num <- matrix(rep(prior,n),n,k,byrow=TRUE)*(matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad) # (n x k)
z.den <- rowSums(z.num) # n-vector
z <- z.num/matrix(rep(z.den,k),ncol=k) # (n x k)
if(is.vector(label))
z[ind.label,] <- class(label,k=k)
numgood <- matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood
numbad <- matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad
v.den <- matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad
v[,1,] <- numgood/v.den
v[,2,] <- numbad/v.den
}
cat("\n")
finalloglik <- loglik[iteration-1]
# ---------------------------------------------------------------------------- #
# The EM-algorithm is finished #
# Check on the EM-monotonicity #
# Plot of the values of the observed-data log-likelihood versus the iterations #
# ---------------------------------------------------------------------------- #
# if(loglikplot){
#
# par(mai=c(0.84,0.8,0.012,0.004))
# par(las = 3)
# par(cex.axis=0.7)
# par(cex.lab=1.2)
# plot(0:(iteration-2),loglik[1:(iteration-1)],type="l",axes = FALSE,xlab="iterations",ylab="log-likelihood",lwd=2)
# axis(1, at = 0:(iteration-2),label = 0:(iteration-2))
# axis(2)
# box(col = "black")
#
# }
# --------------------- #
# Classification Matrix #
# --------------------- #
group <- apply(z,1,which.max)
innergroup <- numeric(n)
for(i in 1:n)
innergroup[i] <- ifelse(v[i,1,group[i]]<v[i,2,group[i]],"bad","*")
detection <- data.frame(group=group,innergroup=innergroup)
# -------------------- #
# Information criteria #
# -------------------- #
if(etafix){
if(!alphafix)
npar <- (k-1) + k + p*k + .ncovpar(modelname=modelname, p=p, G=k)
if(alphafix)
npar <- (k-1) + p*k + .ncovpar(modelname=modelname, p=p, G=k)
}
if(!etafix){
if(!alphafix)
npar <- (k-1) + k + p*k + .ncovpar(modelname=modelname, p=p, G=k) + k
if(alphafix)
npar <- (k-1) + p*k + .ncovpar(modelname=modelname, p=p, G=k) + k
}
AIC <- 2*finalloglik - npar*2
BIC <- 2*finalloglik - npar*log(n)
z.const <- (z<10^(-322))*10^(-322)+(z>10^(-322))*z # vincolo per evitare i NaN nel calcolo di tau*log(tau)
hard.z <- (matrix(rep(apply(z,1,max),k),n,k,byrow=F)==z)*1
if(is.vector(label)){
ECM <- sum(hard.z[-ind.label,]*log(z.const[-ind.label,]))
ICL <- BIC+ECM
}
if(!is.vector(label)){
ECM <- sum(hard.z*log(z.const))
ICL <- BIC+ECM
}
# ---- #
# plot #
# ---- #
# if(plot){
# if(loglikplot)
# x11()
# for(i in 2:p)
# for(j in 1:(i-1)){
# par(mai=c(0.84,0.8,0.012,0.004))
# par(las = 3)
# par(cex.axis=0.7)
# par(cex.lab=1.2)
# plot(X[,c(j,i)],col="white",xlab=paste("X_",j,sep=""),ylab=paste("X_",i,sep=""))
# text(X[,c(j,i)],labels=detection$innergroup,col=group,cex=0.7)
# box(col = "black")
# if((i+j)!=(2*p-1)) x11()
# }
# }
result <- list(
modelname = modelname,
npar = npar,
X = X,
k = k,
p = p,
n = n,
prior = prior,
priorgood = priorgood,
mu = mu,
Sigma = Sigma,
lambda = lambda,
Delta = Delta,
Gamma = Gamma,
eta = eta,
iter.stop = iteration,
z = z,
v = v,
group = group,
detection = detection,
loglik = finalloglik,
AIC = AIC,
BIC = BIC,
ICL = ICL, # alla McNicholas
call = match.call()
)
class(result) <- "MCM"
return(result)
alarm()
}
#####################
## model Selection ##
#####################
MS <- function(
X, # matrix of data
k, # vector of values for k
model=NULL, # models to be considered in model selection
initialization="mclust", # initialization procedure: "random.soft", "random.hard", "manual", or "mclust"
alphacon=TRUE, # if TRUE, alpha is constrained to be >0.5
alphamin=NULL, # minimum proportion of good data in each group
alphafix=FALSE, # if TRUE, the inner weights are estimated
alpha=NULL, # vector of dimension k with proportion of good observations in each group
etacon=TRUE, # if TRUE, eta is constrained to be >1
etafix=FALSE, # if TRUE, eta is not estimated
eta=NULL, # inflation parameters
etamax=200, # maximum value of eta
start.z=NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual"
start.v=NULL, # (n x 2 x k)-array of soft or hard classification in each group: it is used only if initialization="manual"
start=0, # initialization for the package mixture
ind.label=NULL, # indexes of the labelled observations
label=NULL, # groups of the labelled observations
iter.max=1000, # maximum number of iterations in the EM-algorithm
threshold=1.0e-03 # stopping rule in the Aitken rule
)
{
call=match.call()
if (is.data.frame(X))
X <- as.matrix(X)
# if(is.null(G)){
# G=1:3
# }
# if(is.null(modelNames)){
# modelnames=.cwm$modelNames
# }
# if(is.null(method)){
# method="BIC"
# }
n <- length(X)
if(is.null(model))
model <- c("EII","VII","EEI","VEI","EVI","VVI","EEE","VEE","EVE","EEV","VVE","VEV","EVV","VVV")
gridk <- k
numk <- length(gridk)
nummodel <- length(model)
IC <- array(0,c(numk,nummodel,4),dimnames=list(paste(gridk,"groups",sep=" "),paste("model",model,sep=" "),c("loglik","AIC","BIC","ICL")))
cont <- 0
par <- list()
#results <- NULL
#pb <- winProgressBar("Processing Simulations")
#ee <- proc.time()
for(i in 1:numk){
par[[i]] <- list()
#results[[i]] <- NULL
for(j in 1:nummodel){
cat("\n")
cat(paste("Model ",model[j]," with ",gridk[i]," groups",sep=""))
cat("\n")
cont <- cont+1
par[[i]][[j]] <- .MCM(
X=X,
k=gridk[i],
initialization=initialization,
modelname=model[j],
alphacon=alphacon,
alphamin=rep(0.5,gridk[i]),
alphafix=alphafix,
alpha=alpha,
etacon=etacon,
etafix=etafix,
eta=eta,
etamax=etamax,
start.z=start.z,
start.v=start.v,
start=start,
ind.label=ind.label,
label=label,
iter.max=iter.max,
threshold=threshold
#loglikplot=FALSE,
#plot=FALSE
)
#results[[i]][[j]] <- temp
IC[i,j,1] <- par[[i]][[j]]$loglik
IC[i,j,2] <- par[[i]][[j]]$AIC
IC[i,j,3] <- par[[i]][[j]]$BIC
IC[i,j,4] <- par[[i]][[j]]$ICL
#perc <- cont/(numk*nummodel)
#setWinProgressBar(pb,perc)
}
}
#close(pb)
#proc.time()-ee
cat("\n\n")
cat("# ----------------------- #","\n")
cat("# Model Selection Results #","\n")
cat("# ----------------------- #","\n\n")
# --- #
# AIC #
# --- #
BestIndAIC <- which.max(IC[,,2])
BestIndAIC <- arrayInd(BestIndAIC, .dim=c(numk,nummodel))
BestAIC <- IC[BestIndAIC[1],BestIndAIC[2],2]
bestAIC <- par[[BestIndAIC[1]]][[BestIndAIC[2]]]
cat("Best AIC value of",BestAIC,"obtained for k =",gridk[BestIndAIC[1]],"group(s) with model",model[BestIndAIC[2]],"\n\n")
# --- #
# BIC #
# --- #
BestIndBIC <- which.max(IC[,,3])
BestIndBIC <- arrayInd(BestIndBIC, .dim=c(numk,nummodel))
BestBIC <- IC[BestIndBIC[1],BestIndBIC[2],3]
bestBIC <- par[[BestIndBIC[1]]][[BestIndBIC[2]]]
cat("Best BIC value of",BestBIC,"obtained for k =",gridk[BestIndBIC[1]],"group(s) with model",model[BestIndBIC[2]],"\n\n")
# --- #
# ICL #
# --- #
BestIndICL <- which.max(IC[,,4])
BestIndICL <- arrayInd(BestIndICL, .dim=c(numk,nummodel))
BestICL <- IC[BestIndICL[1],BestIndICL[2],4]
bestICL <- par[[BestIndICL[1]]][[BestIndICL[2]]]
cat("Best ICL value of",BestICL,"obtained for k =",gridk[BestIndICL[1]],"group(s) with model",model[BestIndICL[2]],"\n\n")
bestk <- bestmodel <- array(0,c(3),dimnames=list(c("AIC","BIC","ICL")))
bestk <- c(gridk[BestIndAIC[1]],gridk[BestIndBIC[1]],gridk[BestIndICL[1]])
bestmodel <- c(model[BestIndAIC[2]],model[BestIndBIC[2]],model[BestIndICL[2]])
best <- data.frame(k=bestk,model=bestmodel)
return(
structure(
list(
call = call,
best = best,
#IC = IC,
bestAIC = bestAIC,
bestBIC = bestBIC,
bestICL = bestICL
),
class = "MCMMS"
)
)
}
Try the pmcgd package in your browser
Any scripts or data that you put into this service are public.
pmcgd documentation built on May 1, 2019, 7:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions class .cov.wt.punzo .m.step .npar.model .ncovpar .model.type .gpcm .print.gpcm .summary.gpcm .plot.gpcm .igpar .igparv .rgpar .rwgpar .igpark .EM .EMn .EMv .weights .e.step .loglik .getall .MAP .combinewk .msEEE .sumSk.wt .msEEV .getA .getOk .getEkOk .msVEV .msVVV .msEEI .msVEI .msEVI .msVVI .msEII .msVII .msVVE .newD3.MM .newD4.MM .newD .testval .testgrad.D .msEVE .msVVE .msVEE .msEVV .MCM MS
Documented in class MS
#library("mixture")
#library("mnormt")
############################################
## Classification matrix given the groups ##
############################################
class <- function(groups,k){
n <- length(groups)
z <- array(0,c(n,k),dimnames=list(1:n,paste("comp",1:k,sep="")))
for(i in 1:n)
z[i,groups[i]] <- 1
return(z)
}
#################################################################
## Modified Weighted Covariance Matrix for contaminated models ##
#################################################################
.cov.wt.punzo <- function(x,wt,fact){
# fact contains the corrections factors due to the contamination
if (is.data.frame(x))
x <- as.matrix(x)
else if (!is.matrix(x))
stop("'x' must be a matrix or a data frame")
if (!all(is.finite(x)))
stop("'x' must contain finite values only")
p <- ncol(x)
n <- nrow(x)
mu <- array(colSums(wt*fact/sum(wt*fact)*x),c(p),dimnames=list(paste("X.",1:p,sep="")))
cov <- array(crossprod(sqrt(wt*fact/sum(wt))*(x-matrix(rep(mu,n),n,p,byrow=TRUE))),c(p,p),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep="")))
return(
list(
center = mu,
cov = cov
)
)
}
################################
## M-step of the EM algorithm ##
################################
.m.step <- function(data=NULL, covtype=NULL, w=NULL, fact=matrix(1,nrow(w),ncol(w)), v=1, D=NULL, mtol=NULL, mmax=NULL) {
G= ncol(w);
d= ncol(data);
Sk = array(0, c(d,d,G) )
gpar= list()
for (k in 1:G ) {
gpar[[k]] = list()
temp = .cov.wt.punzo(x=data, wt=w[,k], fact=fact[,k])
gpar[[k]]$mu = temp$center
if (!any(is.na( temp$cov))) gpar[[k]]$sigma = temp$cov
Sk[,,k] = temp$cov
}
gpar$pi = apply(w,2,mean)
temp = .model.type(modelname = covtype, Sk=Sk, ng=gpar$pi, D=D, mtol= mtol, mmax= mmax )
gpar$D = temp$D
for (k in 1:G ) {
gpar[[k]]$sigma = temp$sigma[,,k] #+ diag( 1-v , d, d)
gpar[[k]]$invSigma = temp$invSigma[,,k]
gpar[[k]]$logdet = temp$logdet[k]
}
return(gpar)
}
.npar.model <- function(modelname=NULL, p=NULL, G=NULL) {
val = numeric(3)
val[1] = G-1
val[2] = G*p
val[3] = .ncovpar(modelname= modelname, p=p, G=G)
val = sum(val)
return(val)
}
.ncovpar <- function(modelname=NULL, p=NULL, G=NULL) {
if (is.null(p)) stop("p is null")
if (is.null(G)) stop("G is null")
if (is.null(modelname)) stop("modelname is null")
if (modelname == "EII") npar = 1
else if (modelname == "VII") npar = G
else if (modelname == "EEI") npar = p
else if (modelname == "VEI") npar = p + G -1
else if (modelname == "EVI") npar = p*G - G +1
else if (modelname == "VVI") npar = p*G
else if (modelname == "EEE") npar = p*(p+1)/2
else if (modelname == "EEV") npar = G*p*(p+1)/2 - (G-1)*p
else if (modelname == "VEV") npar = G*p*(p+1)/2 - (G-1)*(p-1)
else if (modelname == "VVV") npar = G*p*(p+1)/2
else if (modelname == "EVE") npar = p*(p+1)/2 + (G-1)*(p-1)
else if (modelname == "VVE") npar = p*(p+1)/2 + (G-1)*p
else if (modelname == "VEE") npar = p*(p+1)/2 + (G-1)
else if (modelname == "EVV") npar = G*p*(p+1)/2 - (G-1)
else stop("modelname is not correctly defined")
return(npar)
}
.model.type <- function(modelname=NULL, Sk=NULL, ng=NULL, D=NULL, mtol=1e-10, mmax=10) {
if (is.null(modelname)) stop("modelname is null")
if (modelname == "EII") val = .msEII(Sk=Sk, ng=ng)
else if (modelname == "VII") val = .msVII(Sk=Sk, ng=ng)
else if (modelname == "EEI") val = .msEEI(Sk=Sk, ng=ng)
else if (modelname == "VEI") val = .msVEI(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "EVI") val = .msEVI(Sk=Sk, ng=ng)
else if (modelname == "VVI") val = .msVVI(Sk=Sk, ng=ng)
else if (modelname == "EEE") val = .msEEE(Sk=Sk, ng=ng)
else if (modelname == "EEV") val = .msEEV(Sk=Sk, ng=ng)
else if (modelname == "VEV") val = .msVEV(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "VVV") val = .msVVV(Sk=Sk, ng=ng)
else if (modelname == "EVE") val = .msEVE(Sk=Sk, ng=ng, D0=D, eplison= mtol, max.iter= mmax)
else if (modelname == "VVE") val = .msVVE(Sk=Sk, ng=ng, D0=D, eplison= mtol, max.iter= mmax)
else if (modelname == "VEE") val = .msVEE(Sk=Sk, ng=ng, eplison= mtol, max.iter= mmax)
else if (modelname == "EVV") val = .msEVV(Sk=Sk, ng=ng)
else stop("modelname or covtype is not correctly defined")
if (!is.list(val)) val = list(sigma=val)
return(val)
}
.gpcm <- function(data=NULL, G=1:3, mnames=NULL, start=0, label=NULL, veo=FALSE, nmax=1000, atol=1e-8, mtol=1e-8, mmax=10, pprogress=FALSE, pwarning=TRUE) {
if (is.null(data)) stop('Hey, we need some data, please! data is null')
if (!is.matrix(data)) stop('The data needs to be in matrix form')
if (!is.numeric(data)) stop('The data is required to be numeric')
if (nrow(data) == 1) stop('nrow(data) is equal to 1')
if (ncol(data) == 1) stop('ncol(data) is equal to 1; This function currently only works with multivariate data p > 1')
if (any(is.na(data))) stop('No NAs allowed.')
if (is.null(G)) stop('G is NULL')
G = as.integer(ceiling(G))
if (!is.integer(G)) stop('G is not a integer')
if ( any(G < 1)) stop('G is not a positive integer')
if (is.null(mnames) ) mnames = c("EII", "VII", "EEI", "VEI", "EVI", "VVI", "EEE", "EEV", "VEV", "VVV", "EVE", "VVE", "VEE", "EVV")
bic = array(0, dim= c(length(G), length(mnames), 3), dimnames=list(G, mnames, c('loglik', "npar", "BIC")) )
model = NULL; curBIC = Inf;
for (g in 1:length(G)) {
for (i in 1:length(mnames)) {
if ( pprogress ) print(c(G[g],mnames[i]))
if (veo | .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) < nrow(data)) {
a = try( { .EM(data=data, G=G[g], nmax=nmax, covtype=mnames[i], start=start, label=label, atol= atol, mtol=mtol, mmax=mmax ) }, pwarning)
if ( length(a) > 1) {
bic[g,i,1:2] = c(a$loglik[length(a$loglik)], .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = -2*bic[g,i,1] + bic[g,i,2]*log(nrow(data))
if ( bic[g,i,3] < curBIC) {
model = append(a,list(mtype=mnames[i]))
curBIC = bic[g,i,3]
}
} else {
bic[g,i,1:2] = c(NA, .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = NA
}
} else {
bic[g,i,1:2] = c(NA, .npar.model(modelname=mnames[i], p=ncol(data), G=G[g]) )
bic[g,i,3] = NA
}
}}
if ( is.null(start) ) startobject= "deterministic annealing with default values."
else if ( is.matrix(start)) startobject= "a user specified initialization matrix."
else if ( is.function(start)) startobject= "a user specified initialization function."
else if ( length(start) > 1 ) startobject= "deterministic annealing with user specified values"
else if ( start == 0) startobject= "k-means"
else if ( start > 0) startobject= paste( start, " random initializations", collapse="" )
bicModel = list(G=length(model$gpar$pi), covtype=model$mtype, bic=curBIC )
val = list( start=start, startobject= startobject, gpar=model$gpar, loglik=model$loglik, z=model$z, map=model$MAP, BIC=bic, bicModel= bicModel)
# val = list( model=model, BIC=bic, bicModel= bicModel)
class(val)<-"gpcm"
return(val)
}
.print.gpcm <-function(x, ...){
# a = object$BIC[,,1]
cat("The model choosen by applying the BIC criteria has ", x$bicModel$G, "component(s) and the ", x$bicModel$covtype, "covariance structure\n using ", x$startobject, "\n" )
#endPrint(a)
}
.summary.gpcm <- function(object, ...){
bicl = object$BIC[,,3]
cat("BIC for each model, number of components (rows), and covariance structure (columns).\n")
print.default( bicl)
}
.plot.gpcm <- function(x, ...) {
bicl = x$BIC[,,1]
g = dimnames(bicl)[[1]]
cov = dimnames(bicl)[[2]]
ncov = length(cov)
ng = length(g)
plot( g, bicl[,1], ylim=range(bicl,na.rm=TRUE), xlab="# components", ylab="BIC",type='l')
if (ncov >1) {
for (i in 2:ncov) lines( g, bicl[,i], col=i)
}
legend( "bottomright", legend=cov, col=1:ncov, lty=1)
}
.igpar = function(data=NULL, g=NULL, covtype=NULL, start=NULL, labels=NULL, mtol=NULL, mmax=NULL) {
if (is.null(start)) start = seq( 0.1, 1, length.out=5)
if (is.function(start) ) w = start(data=data, g=g, covtype=covtype)
else if ( !is.null(dim(start)) ) w = start
else if (length(start) > 1) w = .igparv(data=data, g=g, covtype=covtype, vseq0=start , mtol=mtol, mmax = mmax, labels=labels)
else if (start == 0) w = .igpark(data=data, g=g, covtype=covtype)
else if (start > 0) w = .rgpar(data=data, g=g, covtype=covtype, n=ceiling(start), labels=labels, mtol=mtol, mmax = mmax )
else stop(paste('Initialization method ', start, " does not compute!"))
if (!is.matrix(w)) stop("The zij initialization matrix is not a matrix")
if (any(w < 0)) stop("Some of the elements of the zij initialization matrix are less than zero")
# if ( any(apply(w,1,sum) !=1) ) stop("Some of the rows of the zij initialization matrix do not sume to 1")
if (nrow(w) < nrow(data) ) stop("The nrow(zij) of the initialization matrix is less than nrow(data)")
if (nrow(w) > nrow(data) ) stop("The nrow(zij) of the initialization matrix is greater than nrow(data)")
if (ncol(w) < g ) stop("The nrow(zij) of the initialization matrix is less than match g")
if (ncol(w) > g ) stop("The nrow(zij) of the initialization matrix is greater than match g")
w = .combinewk(w, label=labels)
gpar = .m.step(data=data, covtype=covtype, w=w, v=1, mtol=mtol, mmax = mmax)
return(gpar)
}
.igparv <- function(data=NULL, g=NULL,covtype=NULL, vseq0=NULL, labels=NULL, mtol=NULL, mmax = NULL) {
vseq0 = as.numeric(vseq0)
if (is.null(vseq0)) stop('The sequence for deterministic annealing is NULL')
if ( !all( vseq0 <=1 & vseq0 >=0) ) stop('The sequence for deterministic annealing must be between 0 and 1')
w = .rwgpar(data = data, g=g, covtype = covtype)
w = .combinewk(w, label=labels)
gpar = .m.step(data=data, covtype=covtype, w=w, v=1, mtol= mtol, mmax = mmax)
gpar = .EMv(data=data, gpar0 = gpar, G=g, vseq=vseq0, m=1, label=labels, covtype = covtype, mtol= mtol, mmax = mmax )$gpar
w = .weights(data=data, gpar= gpar)
return(w)
}
.rgpar <- function(data=NULL, g=NULL, covtype=NULL, n=1, labels=NULL, mtol=NULL, mmax = NULL) {
w.old = .rwgpar(data= data, g=g,covtype= covtype, labels=labels)
gpar.old = .m.step(data=data, covtype=covtype, w=w.old, v=1, mtol=mtol, mmax = mmax)
loglik.old = .loglik(data=data, gpar=gpar.old)
for (i in 1:n) {
w.new = .rwgpar(data= data, g=g,covtype= covtype, labels=labels)
gpar.new = .m.step(data=data, covtype=covtype, w=w.new, v=1, mtol=mtol, mmax = mmax)
loglik.new = .loglik(data, gpar.new)
if (loglik.new > loglik.old) {
w.old = w.old
gpar.old = gpar.new
loglik.old = loglik.new
}
}
return(w.old)
}
.rwgpar <- function(data=NULL, g=NULL,covtype=NULL, labels=NULL) {
w = matrix(rexp(nrow(data)*g),nrow=nrow(data),ncol=g)
w = matrix(t(apply(w,1, function(z) { z/sum(z)})),nrow=nrow(data),ncol=g)
w = .combinewk(w, label=labels)
return(w)
}
.igpark <- function(data=NULL, g=NULL,covtype=NULL) {
lw = kmeans(data, centers=g, iter.max=10)$cluster
w = .combinewk(matrix(0,nrow=nrow(data),ncol=g), label=lw)
return(w)
}
.EM <- function(data=NULL, gpar0=NULL, G=2, start=1, label=NULL, covtype=NULL, nmax=1000, atol=1e-8, mtol=1e-8, mmax=10 ) {
val = list()
if (is.null(gpar0)) val$gpar = .igpar(data=data, g=G, covtype=covtype, start=start, labels=label, mtol=mtol, mmax=mmax)
else val$gpar = gpar0
val$loglik = numeric(nmax)
val$loglik[1] = .loglik(data=data, gpar=val$gpar)
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[2] = .loglik(data=data, gpar=val$gpar)
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[3] = .loglik(data=data, gpar=val$gpar)
i = 3
while ( ( .getall(val$loglik[(i-2):i]) > atol) & (i < (nmax) ) ) {
i = i+1
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, D=val$gpar$D, mtol=mtol, mmax=mmax)
val$loglik[i] = .loglik(data=data, gpar=val$gpar)
}
val$loglik = val$loglik[1:i]
val$z = .e.step(data=data, gpar=val$gpar, labels=label)
val$MAP = .MAP(data=data, gpar=val$gpar, label=label)
return(val)
}
.EMn <- function(data=NULL, gpar0 = NULL, G=2, n=10, label =NULL, covtype="VVV", mtol=1e-8, mmax=10 ) {
val = list()
if (is.null(gpar0)) val$gpar = .igpar(data=data, g=G, covtype=covtype, mtol=mtol, mmax=mmax)
else val$gpar = gpar0
val$loglik = numeric(n)
for (i in 1:n) {
tempw = .e.step(data=data, gpar=val$gpar, labels=label)
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, mtol=mtol, mmax=mmax)
val$loglik[i] = .loglik(data=data, gpar=val$gpar)
}
# val$label = label
return(val)
}
.EMv <- function(data=NULL, gpar0=NULL, G=3, vseq=c(1,1), m=2, label=NULL, covtype="VVV", mtol=NULL, mmax = NULL ) {
val = list()
if (is.null(gpar0)) val$gpar = .rgpar(data=data, g=G, covtype=covtype, n=1 )
else val$gpar = gpar0
val$loglik = numeric(length(vseq)*m)
count = 1
for (i in 1:length(vseq)) { for (j in 1:m) {
tempw = .e.step(data=data, gpar=val$gpar, labels=label, v=vseq[i])
val$gpar = .m.step(data=data, covtype=covtype, w=tempw, D=val$gpar$D, mtol= mtol, mmax = mmax)
val$loglik[i] = .loglik(data, val$gpar)
count = count + 1
#cat("iteration = ", count-1, "\t v= ", vseq[i], "\t loglik = ", val$loglik[count-1], "\t pi=", val$gpar$pi, "\n")
}}
return(val)
}
.weights <- function(data=NULL, gpar=NULL, v=1) {
d = ncol(data)
G = length(gpar$pi)
if (G > 1) {
zlog = matrix(0, nrow=nrow(data), ncol=length(gpar$pi))
for (k in 1:G ) zlog[,k] = -1/2*mahalanobis(x=data, center=gpar[[k]]$mu, cov=gpar[[k]]$invSigma, inverted=TRUE) -1/2*gpar[[k]]$logdet - d/2 *log(2*pi)
w = t(apply( zlog, 1, function(z,wt,v) {
x= exp( v*(z + log(wt)) )
x=x/sum(x);
return(x) }, wt=gpar$pi,v=v ))
# w = t(apply(zlog, 1, function(z,wt,v) {
# x=exp(v*(z + log(wt)) );
# if (sum(x) == 0) x= rep(1,length(x))
# x = x/sum(x)
# return( x )
# }, wt=gpar$pi,v=v ))
} else w = matrix(1,nrow=nrow(data), ncol=G)
return(w)
}
.e.step <- function(data=NULL, gpar=NULL, labels=NULL, v=1) {
w = .weights(data=data, gpar=gpar,v=v)
if (!is.null(labels)) w = .combinewk(weights=w, label= labels)
return(w)
}
.loglik = function(data, gpar) {
# output is a G x nrow(data) matrix
d = ncol(data)
G = length(gpar$pi)
zlog = matrix(0, nrow=nrow(data), ncol=G)
for (k in 1:G) zlog[,k] = -1/2*mahalanobis(x=data, center=gpar[[k]]$mu, cov=gpar[[k]]$invSigma, inverted=TRUE) -1/2*gpar[[k]]$logdet - d/2*(log(2)+log(pi))
w = apply( exp(zlog),1,function(z,wt) { sum(z*wt) } , wt=gpar$pi)
val = sum(log(w))
if( is.nan(val) ) val =NA
return(val)
}
.getall <- function(loglik) {
lm1 = loglik[3]
lm = loglik[2]
lm_1 = loglik[1]
am = (lm1 - lm)/(lm - lm_1)
lm1.Inf = lm + (lm1 - lm)/(1-am)
val = lm1.Inf - lm
if (is.nan(val) ) val =0
return( abs(val) )
}
.MAP <- function(data, gpar, label=NULL) {
w = .weights(data=data, gpar=gpar, v=1)
if (!is.null(label)) w = .combinewk(weights=w, label= label)
z = apply(w, 1, function(z) { z=(1:length(z))[z==max(z)]; return(z[1]) })
z = as.numeric(z)
return( z)
}
.combinewk <- function(weights=NULL, label=NULL) {
# known is a numeric with
# 0 if unknown group membership
# 1,2,3,.. for label of known group
if (!is.null(label)) { #stop('label is null')
label = as.integer(label)
if (any(!is.integer(label))) stop("Labels are not integers")
if (any(label < 0 )) stop("Labels can only be positive integers")
if (ncol(weights) < max(label) ) stop("Number of groups is less then the number groups given by labels")
if ( sum(label!=0) == nrow(weights) ) {
if (ncol(weights) > max(label) ) stop("Every observations has a label; Cannot fit more groups to the data then given the by the labels.")
}
kw = label !=0
for (j in 1:ncol(weights)) weights[kw,j] = (label == j)[kw]
}
return(weights)
}
.msEEE <- function(Sk=NULL, ng=NULL) {
# Sk is an array of with dim (p x p x G)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = W
logdetW = log(det(W))
invW = solve(W)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = W
val$invSigma[,,g] = invW
val$logdet[g] = logdetW
}
return(val)
}
.sumSk.wt <- function(Sk=NULL, wt=NULL, d=NULL, G=NULL) {
# Sum Sk over the groups with weights wt.
W = matrix(0, nrow=d, ncol=d )
for (g in 1:G) W = W + Sk[,,g]* wt[g]
return(W)
}
.msEEV <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
EWk = Sk
A = matrix(0,d,d)
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk[,,g] = eigen(Wk)$vectors
A = A + t(EWk[,,g]) %*% Wk %*% EWk[,,g]
}
lam = prod(diag(A))^(1/d)
A = A/lam
lam = lam/sum(ng)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam *( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val$invSigma[,,g] = 1/lam *( EWk[,,g] %*% diag(1/diag(A),d) %*% t(EWk[,,g]) )
val$logdet[g] = d*log(lam)
}
return(val)
}
.getA <- function(Ok=NULL, lam=NULL, G=NULL, d=NULL) {
A = matrix(0,d,d)
for (g in 1:G) A = A + Ok[,,g]/lam[g]
A = diag(A)
A = diag( A/prod(A)^(1/d) )
return( A )
}
.getOk <- function(Sk=NULL, ng=NULL, G=NULL) {
Ok = Sk
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk = eigen(Wk)$vectors
Ok[,,g] = t(EWk) %*% Wk %*% EWk
Ok[,,g] = diag(diag(Ok[,,g]))
}
return(Ok)
}
.getEkOk <- function(Sk=NULL, ng=NULL, G=NULL) {
Ok = Sk
EWk = Sk
for (g in 1:G) {
Wk = Sk[,,g]*ng[g]
EWk[,,g] = eigen(Wk)$vectors
Ok[,,g] = t(EWk[,,g]) %*% Wk %*% EWk[,,g]
}
return(list(Ok=Ok,EWk=EWk))
}
.msVEV <- function(Sk=NULL, ng=NULL, eplison=1e-14, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
temp = .getEkOk(Sk=Sk, ng=ng, G=G)
Ok = temp$Ok
EWk = temp$EWk
lam = apply(Ok,3,function(z) { sum(diag(z)) } )/(ng*d)
A = .getA(Ok=Ok, lam=lam, d=d, G=G)
lam = apply(Ok,3,function(z, invA) { sum(diag(z*invA)) }, invA=diag(1/diag(A)) )/(ng*d)
conv = c( d*sum(ng*(1+log(lam))), Inf )
count = 1
while ( diff(conv)/conv[1] > eplison & count < max.iter) {
A = .getA(Ok=Ok, lam=lam, d=d, G=G)
lam = apply(Ok,3,function(z, invA) { sum(diag(z*invA)) }, invA=diag(1/diag(A)) )/(ng*d)
conv = c(d*sum(ng*(1+log(lam))), conv[1] )
count = count +1
}
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = lam[g] * ( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g] * ( EWk[,,g] %*% A %*% t(EWk[,,g]) )
val$invSigma[,,g] = 1/lam[g] * ( EWk[,,g] %*% diag(1/diag(A),d) %*% t(EWk[,,g]) )
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msVVV <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = Sk[,,g]
val$invSigma[,,g] = solve(Sk[,,g])
val$logdet[g] = log(det(Sk[,,g]) )
}
return(val)
}
.msEEI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
B = diag(diag(W))
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = B
val$invSigma[,,g] = diag( 1/diag(B),d )
val$logdet[g] = sum(log( diag(B) ))
}
return(val)
}
.msVEI <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
lam = apply(Sk,3,function(z) { sum(diag(z)) } )/d
# W = sumSk.wt(Sk=Sk, wt=ng*lam, d=d, G=G)
W = .sumSk.wt(Sk=Sk, wt=ng/lam, d=d, G=G)
W = diag(W)
B = diag(W/prod(W)^(1/d))
lam = apply(Sk,3,function(z, invB) { sum(diag(z*invB)) }, invB=diag(1/diag(B) ) )/d
conv = c( d*sum(ng*(1+log(lam))), Inf )
count = 1
while ( abs(diff(conv)) > eplison & count < max.iter) {
# while ( count < max.iter) {
# W = sumSk.wt(Sk=Sk, wt=ng*lam, d=d, G=G)
W = .sumSk.wt(Sk=Sk, wt=ng/lam, d=d, G=G)
W = diag(W)
B = diag(W/prod(W)^(1/d))
lam = apply(Sk,3,function(z, invB) { sum(diag(z*invB)) }, invB=diag(1/diag(B) ) )/d
conv = c(d*sum(ng*(1+log(lam))), conv[1] )
count = count +1
}
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g]*B
val$invSigma[,,g] = diag( 1/diag(B) * 1/lam[g],d )
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msEVI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Bk = matrix(0,nrow=d, ncol=G)
for (g in 1:G) Bk[,g] = diag(Sk[,,g]*ng[g])
lam = apply(Bk, 2, prod)^(1/d)
Bk = sweep(Bk, 2, 1/lam, FUN="*")
lam = sum(lam)/sum(ng)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam * diag(Bk[,g],d)
val$invSigma[,,g] = 1/lam * diag(1/Bk[,g],d)
val$logdet[g] = d*log(lam)
}
return(val)
}
.msVVI <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
# Bk = matrix(0,nrow=d, ncol=G)
# for (g in 1:G) Bk[,g] = diag(Sk[,,g]*ng[g])
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = diag(diag(Sk[,,g]),d)
val$invSigma[,,g] = diag(1/diag(Sk[,,g]),d)
val$logdet[g] = sum( log(diag(Sk[,,g])) )
}
return(val)
}
.msEII <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
W = .sumSk.wt(Sk=Sk, wt=ng, d=d, G=G)/sum(ng)
lam = sum(diag(W))/(sum(ng)*d)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = diag(rep(lam,d), d)
val$invSigma[,,g] = diag(rep(1/lam,d),d)
val$logdet[g] = d*log(lam)
}
return(val)
}
.msVII <- function(Sk=NULL, ng=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
sumdiagSkg = sum(diag(Sk[,,g]))
val$sigma[,,g] = diag(rep(sumdiagSkg/d,d))
val$invSigma[,,g] = diag(rep( d/sumdiagSkg, d))
val$logdet[g] = d* log( sumdiagSkg ) - d*log(d)
}
return(val)
}
.msVVE <- function(Sk=NULL, ng=NULL) {
# Sk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
# lam = numeric(G)
# for (g in 1:G) lam[g] = sum(diag(Sk[,,g]))/d
val = array(0, c(d,d,G))
for (g in 1:G) val[,,g] = diag(rep(sum(diag(Sk[,,g]))/d,d))
return(val)
}
#####################
### new.cov models
#####################
.newD3.MM <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
z = matrix(0,d,d)
lambda =0
for (g in 1:G) {
lambdak = max(eigen(Wk[,,g])$values)
z = z + diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g] -lambdak *(diag(1/Ak[,g])%*% t(D) )
}
z1 = svd(z)
Xk1 = (z1$v) %*% t(z1$u)
return( Xk1 )
}
.newD4.MM <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
z = matrix(0,d,d)
lambda =0
for (g in 1:G) {
lambdak = max(1/Ak[,g])
#z = z + diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g] -lambdak *(diag(1/Ak[,g])%*% t(D) )
z = z + Wk[,,g] %*% (D) %*% diag(1/Ak[,g]) -lambdak *(Wk[,,g]%*% (D) )
}
z1 = svd(z)
# Xk1 = (z1$v) %*% t(z1$u)
# return( t(Xk1) )
# OR
Xk1 = (z1$v) %*% t(z1$u)
return( t(Xk1) )
}
.newD <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL, tmax = 100) {
D6 = D
D6 = .newD3.MM(D=D6, d=d, G=G, Wk=Wk, Ak=Ak, tmax = 100)
D6 = .newD4.MM(D=D6, d=d, G=G, Wk=Wk, Ak=Ak, tmax = 100)
return(D6)
}
.testval <- function(Wk=NULL, Ak=NULL, D=NULL, G=NULL) {
z = numeric(G)
# for (g in 1:G) z[g] = sum(diag( D %*% diag(1/Ak[,g]) %*% t(D) %*% Wk[,,g]))
for (g in 1:G) z[g] = sum(diag( t(D) %*% Wk[,,g] %*% D %*% diag(1/Ak[,g]) ))
return(sum(z))
}
.testgrad.D <- function(D=NULL, d=NULL, G=NULL, Wk=NULL, Ak=NULL) {
z = matrix(0,d,d)
for (g in 1:G) z = z + Wk[,,g] %*% D %*% diag(1/Ak[,g])
return(2*z)
}
.msEVE <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 10, D0=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1]; G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
# D = diag(rep(1,d))
D = D0
if (is.null(D)) D = diag(rep(1,d))
# D = t(eigen(W)$vectors)
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% (D) ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(0, 1, testval(Wk=Wk,Ak=Ak,D=D,G=G)))
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
#print(c(0, 2, testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c( .testval(Wk=Wk,Ak=Ak,D=D,G=G), Inf )
count = 1
while ( diff(conv)/abs(conv[1]) > eplison & count < max.iter) {
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% (D) ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(count, 0, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c(.testval(Wk=Wk,Ak=Ak,D=D,G=G), conv[1] )
count = count +1
}
lam = 0
for (g in 1:G) lam = lam + sum(diag( D %*% diag(1/Ak[,g])%*% t(D) %*% Wk[,,g] ))
lam = lam/(sum(ng)*d)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G), D=D )
for (g in 1:G) {
val$sigma[,,g] = D %*% diag(lam*Ak[,g])%*% t(D)
val$invSigma[,,g] = D %*% diag(1/lam*1/Ak[,g])%*% t(D)
val$logdet[g] = d*log( lam )
}
return(val)
}
.msVVE <- function(Sk=NULL, ng=NULL, eplison=1e-20, max.iter= 10, D0=NULL) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1]; G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
D = D0
if (is.null(D)) D = diag(rep(1,d))
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% D ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
#print(c(0, 1, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
#print(c(0, 2, .testval(Wk=Wk,Ak=Ak,D=D,G=G)))
conv = c( .testval(Wk=Wk,Ak=Ak,D=D,G=G), Inf )
count = 1
while ( diff(conv)/abs(conv[1]) > eplison & count < max.iter) {
Ak = apply(Wk,3, function(z,D) { diag( t(D) %*% z %*% D ) }, D=D )
Ak = apply(Ak,2,function(z) { z/prod(z)^(1/length(z))})
D = .newD(D=D, Wk=Wk, Ak=Ak, G=G, d=d)
conv = c(.testval(Wk=Wk,Ak=Ak,D=D,G=G), conv[1] )
count = count +1
}
#print(count)
lam = numeric(G)
for (g in 1:G) lam[g] =sum(diag( D %*% diag(1/Ak[,g])%*% t(D) %*% Sk[,,g] ))/d
#print(apply(Ak,2,prod))
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = (D %*% diag(lam[g]*Ak[,g])%*% t(D))
val$invSigma[,,g] = (D %*% diag(1/lam[g]*1/Ak[,g])%*% t(D))
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msVEE <- function(Sk=NULL, ng=NULL, eplison=1e-14, max.iter=100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Wk = Sk
W = matrix(0,d,d)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
W = W + Wk[,,g]
}
C = W/det(W)^(1/d)
invC = solve(C)
lam = apply(Sk,3,function(z, invC) { sum(diag(z %*% invC)) }, invC=invC)/d
val1 = sum(apply(Wk,3,function(z, invC) { sum(diag(z*invC)) }, invC= invC)/lam) +d*sum(ng*lam)
conv = c(val1, Inf )
count = 1
while ( diff(conv)/conv[1] > eplison & count < max.iter) {
#for (i in 1:max.iter) {
C = .sumSk.wt(Sk=Wk, wt=1/lam, d=d,G=G)
C = C/det(C)^(1/d)
invC = solve(C)
lam = apply(Sk,3,function(z, invC) { sum(diag(z %*% invC)) }, invC=invC)/d
val1 = sum(apply(Wk,3,function(z, invC) { sum(diag(z*invC)) }, invC= invC )/lam) +d*sum(ng*lam)
conv = c(val1, conv[1] )
count = count +1
}
#print(c(count,det(C), lam) )
invC = solve(C)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam[g]*C
val$invSigma[,,g] = 1/lam[g]*invC
val$logdet[g] = d*log(lam[g])
}
return(val)
}
.msEVV <- function(Sk=NULL, ng=NULL, eplison=1e-12, max.iter= 100) {
# Wk is a list of length G with matrix (p x p)
# ng is a vector of length G of weights of Wk
d = dim(Sk)[1];
G = dim(Sk)[3];
Wk = Sk
Ck = Sk
lam = numeric(G)
for (g in 1:G) {
Wk[,,g] = Sk[,,g]*ng[g]
lam[g] = det(Wk[,,g])^(1/d)
}
Ck = sweep(Wk, 3, 1/lam, FUN="*")
lam = sum(lam)
val = list(sigma=array(0, c(d,d,G)), invSigma=array(0, c(d,d,G)), logdet=numeric(G) )
for (g in 1:G) {
val$sigma[,,g] = lam * Ck[,,g]
val$invSigma[,,g] = 1/lam * solve(Ck[,,g])
val$logdet[g] = d* log(lam )
}
return(val)
}
####################################
## Mixture of Contaminated Models ##
####################################
.MCM <- function(
X, # matrix of data
k, # number of groups
initialization="mclust", # initialization procedure: "random.soft", "random.hard", "manual", or "mclust" (which is nested)
modelname="VVV", # one of the 14 models of Celeaux & Govaert (1995)
alphacon=TRUE, # if TRUE, alpha is constrained to be >0.5
alphamin=rep(0.5,k), # minimum proportion of good data in each group
alphafix=FALSE, # if TRUE, the inner weights are not estimated
alpha=NULL, # vector of dimension k with proportion of good observations in each group
etacon=TRUE, # if TRUE, eta is constrained to be >1
etafix=FALSE, # if TRUE, eta is not estimated
eta=NULL, # inflation parameters
etamax=200, # maximum value of eta
start.z=NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual"
start.v=NULL, # (n x 2 x k)-array of soft or hard classification in each group: it is used only if initialization="manual"
start=0, # initialization for the package mixture
ind.label=NULL, # indexes of the labelled observations
label=NULL, # groups of the labelled observations
iter.max=1000, # maximum number of iterations in the EM-algorithm
threshold=1.0e-04 # stopping rule in the Aitken rule
#loglikplot=TRUE, # if TRUE, the log-likelihood values against the iterations are plotted
#plot=TRUE # if TRUE, the plot of the classified data is given
)
{
#call = match.call()
if (is.data.frame(X))
X <- as.matrix(X)
n <- nrow(X) # sample size
p <- ncol(X) # number of variables
if (is.null(X)) stop('Hey, we need some data, please! X is null')
if (!is.matrix(X)) stop('X needs to be in matrix form')
if (!is.numeric(X)) stop('X is required to be numeric')
if (n == 1) stop('nrow(X) is equal to 1')
if (p == 1) stop('ncol(X) is equal to 1; This function currently only works with multivariate data p > 1')
if (any(is.na(X))) stop('No NAs allowed.')
if (is.null(k)) stop('k is NULL')
k <- as.integer(ceiling(k))
if (!is.integer(k)) stop('k is not a integer')
if (any(k < 1)) stop('k is not a positive integer')
# for model-based classification
lab <- NULL
if(is.vector(label)){
nlab <- length(label)
nunlab <- n-nlab
lab <- numeric(n)
lab[ind.label] <- label
}
prior <- numeric(k) # proportion of each group
priorgood <- numeric(k) # proportion of good observations in each group
v <- array(0.5,c(n,2,k),dimnames=list(1:n,c("good","bad"),paste("group ",1:k,sep="")))
mu <- array(0,c(p,k),dimnames=list(paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
Sigma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
invSigma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
lambda <- array(0,c(k),dimnames = list(paste("group ",1:k,sep="")))
Delta <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
Gamma <- array(0,c(p,p,k),dimnames=list(paste("X.",1:p,sep=""),paste("X.",1:p,sep=""),paste("group ",1:k,sep="")))
if(!etafix)
eta <- rep(2,k) # inflation parameters
correction <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep=""))) # factor which differentiates this model by mclust
PXgood <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep="")))
PXbad <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep="")))
# ------------------------- #
# posteriors initialization #
# ------------------------- #
if(initialization=="random.soft"){
z <- array(runif(n*k),c(n,k)) # soft posterior probabilities (no-normalized) (n x k)
z <- z/rowSums(z) # soft posterior probabilities (n x k)
for(j in 1:k){
v[,2,j] <- runif(n,0,0.2)
v[,1,j] <- 1-v[,2,j]
}
}
if(initialization=="random.hard"){
z <- t(rmultinom(n, size = 1, prob=rep(1/k,k))) # hard posterior probabilities (n x k)
for(j in 1:k)
v[,,j] <- t(rmultinom(n, size = 1, prob=c(0.9,0.1)))
}
if(initialization=="manual"){
z <- start.z # soft or hard classification
v <- start.v # soft or hard classification of good and bad observations
}
if(initialization=="mclust"){
mclustfit <- .gpcm(data=X, G=k, mnames=modelname, start=start, label=lab, pprogress=FALSE)
for(j in 1:k){
mu[,j] <- mclustfit$gpar[[j]]$mu
Sigma[,,j] <- mclustfit$gpar[[j]]$sigma
invSigma[,,j] <- mclustfit$gpar[[j]]$invSigma
#temp <- eigen(Sigma[,,j])
#Gamma[,,j] <- temp$vectors
#Delta[,,j] <- diag(temp$values)/prod(temp$values)^(1/p)
#lambda[j] <- prod(temp$values)^(1/p)
}
z <- mclustfit$z
group <- apply(z,1,which.max)
for(j in 1:k){
Probj <- (2*pi)^(-p/2)*det(Sigma[,,j])^(-1/2)*exp(-1/2*mahalanobis(x=X[which(group==j),], center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
#Probj <- dmnorm(x=X[which(group==j),], mean = mu[,j], varcov=Sigma[,,j])
wgood <- Probj/max(Probj)
v[which(group==j),1,j] <- (wgood)^0.2
v[which(group==j),2,j] <- 1-(wgood)^0.2
}
}
if(is.vector(label))
z[ind.label,] <- class(label,k=k)
# ------------ #
# EM algorithm #
# ------------ #
# Preliminary definition of convergence criterions
check <- 0
iteration <- 1
loglik <- NULL
aloglik <- NULL
aloglik <- c(0,0)
a <- NULL
a <- c(0,0)
while(check<1){
# ++++++ #
# M-step #
# ++++++ #
# ------- #
# Weights #
# ------- #
prior <- colMeans(z)
zv <- z*v[,1,]
zvcompl <- z*v[,2,]
if(!alphafix){
if(alphacon){
g <- function(alpha,j,z,v)
sum(z[,j]*(v[,1,j]*log(alpha)+v[,2,j]*log(1-alpha)))
for(j in 1:k)
priorgood[j] <- optimize(g,c(alphamin[j],1),maximum = TRUE,j=j,z=z,v=v)$maximum
}
if(!alphacon)
priorgood <- colSums(zv)/colSums(z)
}
if(alphafix)
priorgood <- alpha
priorbad <- 1-priorgood
# ---------- #
# mu & Sigma #
# ---------- #
correction <- v[,1,]+v[,2,]*matrix(rep(1/eta),n,k,byrow=TRUE)
fitM <- .m.step(data=X, covtype=modelname, w=z, fact=correction, v=1, mtol=1e-10, mmax=10)
for(j in 1:k){
mu[,j] <- fitM[[j]]$mu
Sigma[,,j] <- fitM[[j]]$sigma
invSigma[,,j] <- fitM[[j]]$invSigma
temp <- eigen(Sigma[,,j])
Gamma[,,j] <- temp$vectors
Delta[,,j] <- diag(temp$values)/prod(temp$values)^(1/p)
lambda[j] <- prod(temp$values)^(1/p)
}
# ------------------- #
# Inflation parameter #
# ------------------- #
if(!etafix){
if(etacon){
f <- function(eta,j,p,zvcompl,X,mu,invSigma)
sum(-(p/2)*zvcompl[,j]*log(eta)-(1/2)*zvcompl[,j]*(1/eta)*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
for(j in 1:k)
eta[j] <- optimize(f,c(1,etamax),maximum = TRUE,j=j,p=p,zvcompl=zvcompl,X=X,mu=mu,invSigma=invSigma)$maximum
}
if(!etacon){
for(j in 1:k)
eta[j] <- sum(zvcompl[,j]*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))/(p*sum(zvcompl[,j]))
}
}
# ------- #
# density #
# ------- #
zerostar <- 1e-323 # to avoid zero probabilities
for(j in 1:k){
PXgood[,j] <- (2*pi)^(-p/2)*(det(Sigma[,,j]))^(-1/2)*exp(-1/2*mahalanobis(x=X, center=mu[,j], cov=invSigma[,,j], inverted=TRUE))
PXbad[,j] <- (2*pi)^(-p/2)*(det(eta[j]*Sigma[,,j]))^(-1/2)*exp(-1/2*mahalanobis(x=X, center=mu[,j], cov=1/eta[j]*invSigma[,,j], inverted=TRUE))
PXgood[,j] <- (PXgood[,j]<zerostar)*zerostar+(PXgood[,j]>=zerostar)*PXgood[,j]
PXbad[,j] <- (PXbad[,j]<zerostar)*zerostar+(PXbad[,j]>=zerostar)*PXbad[,j]
}
# ------------------------------------- #
# Global - Observed-data log-likelihood #
# ------------------------------------- #
# model-based clustering
if(!is.vector(label))
llvalues <- sum(log( rowSums(matrix(rep(prior,n),n,k,byrow=TRUE)*(matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad))))
# model-based classification
if(is.vector(label)){
llvalueslab <- z[ind.label,]*(log(matrix(rep(prior,nlab),nlab,k,byrow=T))+log( matrix(rep(priorgood,nlab),nlab,k,byrow=TRUE)*PXgood[ind.label,]+matrix(rep(priorbad,nlab),nlab,k,byrow=TRUE)*PXbad[ind.label,] ))
llvaluesunlab <- log(rowSums(matrix(rep(prior,nunlab),nunlab,k,byrow=TRUE)*(matrix(rep(priorgood,nunlab),nunlab,k,byrow=TRUE)*PXgood[-ind.label,]+matrix(rep(priorbad,nunlab),nunlab,k,byrow=TRUE)*PXbad[-ind.label,])))
llvalues <- sum(llvalueslab)+sum(llvaluesunlab)
}
loglik[iteration] <- llvalues
# ----------------------------------------------- #
# Aitkane's Acceleration-Based Stopping Criterion #
# ----------------------------------------------- #
if(iteration>2 & k > 1){
if(abs(loglik[iteration-1]-loglik[iteration-2])>0){
a[iteration-1] <- (loglik[iteration]-loglik[iteration-1])/(loglik[iteration-1]-loglik[iteration-2])
aloglik[iteration] <- loglik[iteration-1]+(1/(1-a[iteration-1])*(loglik[iteration]-loglik[iteration-1]))
if(abs(aloglik[iteration]-loglik[iteration])<threshold)
check <- 1
}
#if(abs(loglik[iteration-1]-loglik[iteration-2])==0)
else
check <- 1
}
if(iteration==iter.max | k==1) check <- 1
cat("*")
#cat("Iteration",iteration,"\n")
iteration <- iteration + 1
# ++++++ #
# E-Step #
# ++++++ #
z.num <- matrix(rep(prior,n),n,k,byrow=TRUE)*(matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad) # (n x k)
z.den <- rowSums(z.num) # n-vector
z <- z.num/matrix(rep(z.den,k),ncol=k) # (n x k)
if(is.vector(label))
z[ind.label,] <- class(label,k=k)
numgood <- matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood
numbad <- matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad
v.den <- matrix(rep(priorgood,n),n,k,byrow=TRUE)*PXgood+matrix(rep(priorbad,n),n,k,byrow=TRUE)*PXbad
v[,1,] <- numgood/v.den
v[,2,] <- numbad/v.den
}
cat("\n")
finalloglik <- loglik[iteration-1]
# ---------------------------------------------------------------------------- #
# The EM-algorithm is finished #
# Check on the EM-monotonicity #
# Plot of the values of the observed-data log-likelihood versus the iterations #
# ---------------------------------------------------------------------------- #
# if(loglikplot){
#
# par(mai=c(0.84,0.8,0.012,0.004))
# par(las = 3)
# par(cex.axis=0.7)
# par(cex.lab=1.2)
# plot(0:(iteration-2),loglik[1:(iteration-1)],type="l",axes = FALSE,xlab="iterations",ylab="log-likelihood",lwd=2)
# axis(1, at = 0:(iteration-2),label = 0:(iteration-2))
# axis(2)
# box(col = "black")
#
# }
# --------------------- #
# Classification Matrix #
# --------------------- #
group <- apply(z,1,which.max)
innergroup <- numeric(n)
for(i in 1:n)
innergroup[i] <- ifelse(v[i,1,group[i]]<v[i,2,group[i]],"bad","*")
detection <- data.frame(group=group,innergroup=innergroup)
# -------------------- #
# Information criteria #
# -------------------- #
if(etafix){
if(!alphafix)
npar <- (k-1) + k + p*k + .ncovpar(modelname=modelname, p=p, G=k)
if(alphafix)
npar <- (k-1) + p*k + .ncovpar(modelname=modelname, p=p, G=k)
}
if(!etafix){
if(!alphafix)
npar <- (k-1) + k + p*k + .ncovpar(modelname=modelname, p=p, G=k) + k
if(alphafix)
npar <- (k-1) + p*k + .ncovpar(modelname=modelname, p=p, G=k) + k
}
AIC <- 2*finalloglik - npar*2
BIC <- 2*finalloglik - npar*log(n)
z.const <- (z<10^(-322))*10^(-322)+(z>10^(-322))*z # vincolo per evitare i NaN nel calcolo di tau*log(tau)
hard.z <- (matrix(rep(apply(z,1,max),k),n,k,byrow=F)==z)*1
if(is.vector(label)){
ECM <- sum(hard.z[-ind.label,]*log(z.const[-ind.label,]))
ICL <- BIC+ECM
}
if(!is.vector(label)){
ECM <- sum(hard.z*log(z.const))
ICL <- BIC+ECM
}
# ---- #
# plot #
# ---- #
# if(plot){
# if(loglikplot)
# x11()
# for(i in 2:p)
# for(j in 1:(i-1)){
# par(mai=c(0.84,0.8,0.012,0.004))
# par(las = 3)
# par(cex.axis=0.7)
# par(cex.lab=1.2)
# plot(X[,c(j,i)],col="white",xlab=paste("X_",j,sep=""),ylab=paste("X_",i,sep=""))
# text(X[,c(j,i)],labels=detection$innergroup,col=group,cex=0.7)
# box(col = "black")
# if((i+j)!=(2*p-1)) x11()
# }
# }
result <- list(
modelname = modelname,
npar = npar,
X = X,
k = k,
p = p,
n = n,
prior = prior,
priorgood = priorgood,
mu = mu,
Sigma = Sigma,
lambda = lambda,
Delta = Delta,
Gamma = Gamma,
eta = eta,
iter.stop = iteration,
z = z,
v = v,
group = group,
detection = detection,
loglik = finalloglik,
AIC = AIC,
BIC = BIC,
ICL = ICL, # alla McNicholas
call = match.call()
)
class(result) <- "MCM"
return(result)
alarm()
}
#####################
## model Selection ##
#####################
MS <- function(
X, # matrix of data
k, # vector of values for k
model=NULL, # models to be considered in model selection
initialization="mclust", # initialization procedure: "random.soft", "random.hard", "manual", or "mclust"
alphacon=TRUE, # if TRUE, alpha is constrained to be >0.5
alphamin=NULL, # minimum proportion of good data in each group
alphafix=FALSE, # if TRUE, the inner weights are estimated
alpha=NULL, # vector of dimension k with proportion of good observations in each group
etacon=TRUE, # if TRUE, eta is constrained to be >1
etafix=FALSE, # if TRUE, eta is not estimated
eta=NULL, # inflation parameters
etamax=200, # maximum value of eta
start.z=NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual"
start.v=NULL, # (n x 2 x k)-array of soft or hard classification in each group: it is used only if initialization="manual"
start=0, # initialization for the package mixture
ind.label=NULL, # indexes of the labelled observations
label=NULL, # groups of the labelled observations
iter.max=1000, # maximum number of iterations in the EM-algorithm
threshold=1.0e-03 # stopping rule in the Aitken rule
)
{
call=match.call()
if (is.data.frame(X))
X <- as.matrix(X)
# if(is.null(G)){
# G=1:3
# }
# if(is.null(modelNames)){
# modelnames=.cwm$modelNames
# }
# if(is.null(method)){
# method="BIC"
# }
n <- length(X)
if(is.null(model))
model <- c("EII","VII","EEI","VEI","EVI","VVI","EEE","VEE","EVE","EEV","VVE","VEV","EVV","VVV")
gridk <- k
numk <- length(gridk)
nummodel <- length(model)
IC <- array(0,c(numk,nummodel,4),dimnames=list(paste(gridk,"groups",sep=" "),paste("model",model,sep=" "),c("loglik","AIC","BIC","ICL")))
cont <- 0
par <- list()
#results <- NULL
#pb <- winProgressBar("Processing Simulations")
#ee <- proc.time()
for(i in 1:numk){
par[[i]] <- list()
#results[[i]] <- NULL
for(j in 1:nummodel){
cat("\n")
cat(paste("Model ",model[j]," with ",gridk[i]," groups",sep=""))
cat("\n")
cont <- cont+1
par[[i]][[j]] <- .MCM(
X=X,
k=gridk[i],
initialization=initialization,
modelname=model[j],
alphacon=alphacon,
alphamin=rep(0.5,gridk[i]),
alphafix=alphafix,
alpha=alpha,
etacon=etacon,
etafix=etafix,
eta=eta,
etamax=etamax,
start.z=start.z,
start.v=start.v,
start=start,
ind.label=ind.label,
label=label,
iter.max=iter.max,
threshold=threshold
#loglikplot=FALSE,
#plot=FALSE
)
#results[[i]][[j]] <- temp
IC[i,j,1] <- par[[i]][[j]]$loglik
IC[i,j,2] <- par[[i]][[j]]$AIC
IC[i,j,3] <- par[[i]][[j]]$BIC
IC[i,j,4] <- par[[i]][[j]]$ICL
#perc <- cont/(numk*nummodel)
#setWinProgressBar(pb,perc)
}
}
#close(pb)
#proc.time()-ee
cat("\n\n")
cat("# ----------------------- #","\n")
cat("# Model Selection Results #","\n")
cat("# ----------------------- #","\n\n")
# --- #
# AIC #
# --- #
BestIndAIC <- which.max(IC[,,2])
BestIndAIC <- arrayInd(BestIndAIC, .dim=c(numk,nummodel))
BestAIC <- IC[BestIndAIC[1],BestIndAIC[2],2]
bestAIC <- par[[BestIndAIC[1]]][[BestIndAIC[2]]]
cat("Best AIC value of",BestAIC,"obtained for k =",gridk[BestIndAIC[1]],"group(s) with model",model[BestIndAIC[2]],"\n\n")
# --- #
# BIC #
# --- #
BestIndBIC <- which.max(IC[,,3])
BestIndBIC <- arrayInd(BestIndBIC, .dim=c(numk,nummodel))
BestBIC <- IC[BestIndBIC[1],BestIndBIC[2],3]
bestBIC <- par[[BestIndBIC[1]]][[BestIndBIC[2]]]
cat("Best BIC value of",BestBIC,"obtained for k =",gridk[BestIndBIC[1]],"group(s) with model",model[BestIndBIC[2]],"\n\n")
# --- #
# ICL #
# --- #
BestIndICL <- which.max(IC[,,4])
BestIndICL <- arrayInd(BestIndICL, .dim=c(numk,nummodel))
BestICL <- IC[BestIndICL[1],BestIndICL[2],4]
bestICL <- par[[BestIndICL[1]]][[BestIndICL[2]]]
cat("Best ICL value of",BestICL,"obtained for k =",gridk[BestIndICL[1]],"group(s) with model",model[BestIndICL[2]],"\n\n")
bestk <- bestmodel <- array(0,c(3),dimnames=list(c("AIC","BIC","ICL")))
bestk <- c(gridk[BestIndAIC[1]],gridk[BestIndBIC[1]],gridk[BestIndICL[1]])
bestmodel <- c(model[BestIndAIC[2]],model[BestIndBIC[2]],model[BestIndICL[2]])
best <- data.frame(k=bestk,model=bestmodel)
return(
structure(
list(
call = call,
best = best,
#IC = IC,
bestAIC = bestAIC,
bestBIC = bestBIC,
bestICL = bestICL
),
class = "MCMMS"
)
)
}
Try the pmcgd package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.