Nothing
##################################
# DEPMIX INTERNAL FUNCTIONS, NOT TO BE CALLED BY USER
#
#
# AUXILIARIES FOR FITDMM
#
checkSetRecode <- function(dat,dmm,tdcov=0,printlevel=1) {
## create the model to be fitted, changed into class mgd
if(class(dmm)[1]=="dmm") {
xgmod <- mixdmm(dmm=list(dmm))
xgmod <- mgdmm(dmm=xgmod)
}
if(class(dmm)[1]=="mixdmm") xgmod <- mgdmm(dmm=dmm)
if(class(dmm)[1]=="mgd") xgmod <- dmm
## some correctness checks, far from complete
if(xgmod$ng==1) {
if(!(class(dat)=="md")) stop("Data should have class md")
dat <- list(dat)
}
if(length(dat)!=xgmod$ng) stop("Number of data sets does not match nr of groups.")
if(!(length(xgmod$nstates)==xgmod$nrcomp)) stop("Length of nstates and nrcomp do not match.")
if(!(length(xgmod$itemtypes)==(length(itemtypes(dat[[1]]))-ncov(dat[[1]])))) stop("Model incompatible with data, different numbers of items.")
ncov=ncov(dat[[1]])
if(tdcov==1 && ncov==0) stop("There are no covariates in the data so these cannot be used.")
if(tdcov) if(xgmod$td==FALSE) stop("Model does not specify covariates.")
if(printlevel>19) print("All checks okay.")
gc()
## split data and covariates, data may have covariates, even when they're not fitted,
## so the data allways has to be split
if(ncov>0) {
nrit=dim(dat[[1]])[2]-ncov
itt=rep(1,ncov)
dcov=list()
datsplit=list()
for(i in 1:xgmod$ng) {
dcov[[i]]=markovdata(dat=dat[[i]][,(nrit+1):(dim(dat[[i]])[2]),drop=FALSE],itemtypes=itt,ntimes=ntimes(dat[[i]]))
datsplit[[i]]=markovdata(dat=dat[[i]][,1:nrit,drop=FALSE],itemtypes=itemtypes(dat[[i]])[1:nrit],ntimes=ntimes(dat[[i]]),replicates=replicates(dat[[i]]))
}
dat=datsplit
if(printlevel>19) print("Data split okay.")
}
## check missing data codes
xm=-9999
for(i in 1:xgmod$ng) {
dat[[i]]=replace(dat[[i]],which(is.na(dat[[i]])),xm)
if(ncov>0) datsplit[[i]]=replace(datsplit[[i]],which(is.na(datsplit[[i]])),xm)
}
if(printlevel>19) print("Missing data codes set.")
## send td model to C routines
if(tdcov) {
nitems=length(xgmod$itemtypes)
z <- .C("covSetUp",
as.integer(xgmod$ng),
as.integer(xgmod$nrcomp),
as.integer(xgmod$nstates),
as.integer(nitems),
as.integer(xgmod$itemtypes),
as.double(xgmod$pars[(xgmod$npars+1):(2*xgmod$npars)]),
as.integer(xm),
as.integer(printlevel),
PACKAGE="depmix")
## ... and set covariates in C
mcov <- .C("multiCovSetUp",
as.integer(length(dcov)),
as.integer(printlevel),
PACKAGE="depmix")
datset <- list()
for(i in 1:xgmod$ng) {
# catit=0
# dc=numeric(0)
datset[[i]] <- .C("ngCovSetUp",
as.integer(i),
as.double(t(dcov[[i]])),
as.integer(dim(dcov[[i]])[2]),
as.integer(itemtypes(dcov[[i]])),
as.integer(length(ntimes(dcov[[i]]))),
as.integer(ntimes(dcov[[i]])),
as.double(replicates(dcov[[i]])),
as.integer(xm),
as.integer(printlevel),
PACKAGE="depmix")
}
if(printlevel>19) print("Td model set in C")
}
nitems=length(xgmod$itemtypes)
pars=xgmod$pars[1:xgmod$npars] # only the normal model pars, tdpars are already set
if(printlevel>19) print("Preparing to send model to C routines.")
## send model to C routines
z <- .C("mixModelSetUp",
as.integer(xgmod$ng),
as.integer(xgmod$nrcomp),
as.integer(xgmod$nstates),
as.integer(nitems),
as.integer(xgmod$itemtypes),
as.double(pars),
as.integer(xm),
as.integer(printlevel),
PACKAGE="depmix")
if(printlevel>19) print("Model set up finished in C.")
## send data to C routines
mdat <- .C("multiDataSetUp",
as.integer(length(dat)),
as.integer(printlevel),
PACKAGE="depmix")
datset <- list()
for(i in 1:xgmod$ng) {
attr(dat[[i]],"itemtypes")=replace(itemtypes(dat[[i]]),which(itemtypes(dat[[i]])!="categorical"),1)
itt=rep(0,dim(dat[[i]])[2])
## recode itemtypes
for(it in 1:dim(dat[[i]])[2]) {
if(itemtypes(dat[[i]])[it]=="categorical") itt[it]=length(setdiff(unique(dat[[i]][,it]),xm))
else itt[it]=1
}
attr(dat[[i]],"itemtypes")=itt
# recode the data if necessary
catit=0
dc=numeric(0)
for(it in 1:dim(dat[[i]])[2]) {
if(itemtypes(dat[[i]])[it]>1) {
catit=catit+1
datacats=sort(setdiff(unique(dat[[i]][,it]),xm))
if(all.equal(datacats,(1:itt[it]))==TRUE) dc[catit]=1
else dc[catit]=0
}
}
if(all(as.logical(dc))) recdat=dat[[i]]
else recdat=recode(dat[[i]],xm)
datset[[i]] <- .C("ngDataSetUp",
as.integer(i),
as.double(t(recdat)),
as.integer(dim(dat[[i]])[2]),
as.integer(itemtypes(dat[[i]])),
as.integer(length(ntimes(dat[[i]]))),
as.integer(ntimes(dat[[i]])),
as.double(replicates(dat[[i]])),
as.integer(xm),
as.integer(printlevel),
PACKAGE="depmix")
}
if(printlevel>19) print("Data set up finished in C\n")
return(xgmod)
}
recode <- function(dat,xm) {
recdat=dat
nit=length(itemtypes(dat))
for(it in 1:nit) {
if(itemtypes(dat)[it]>1) {
original=sort(setdiff(unique(dat[,it]),xm))
normvals=1:itemtypes(dat)[it]
for(i in 1:itemtypes(dat)[it]) {
recdat[which(dat[,it]==original[i]),it]=normvals[i]
}
}
}
recdat
}
#
# AUXILIARIES FOR DMM
#
# function fblo gives lower bounds of parameters for each density
fblo<-function(x,i,bigB) {
if(x > 1) return(rep(0,x)) # multinomial
if(x== 1) return(c(-bigB,-bigB)) # normal (mu - mean, sigma - stddev)
}
# function fbuo gives upper bounds of parameters for each density
fbuo<-function(x,i,bigB) {
ubshift=+bigB
if(x > 1) return(rep(1,x)) # multinomial
if(x==1) return(c(+bigB,+bigB)) # normal (mu - mean, sigma - stddev)
}
# function ppar generates random start values of parameters for each density and pastes stvals together
ppar<-function(x,z) {
if(x > 1) {y<-runif(x);y<-y/sum(y)} # multinomial
if(x==1) y<-c(rnorm(1), rnorm(1,1,0.1)) # normal (mu - mean, sigma - stddev)
return(c(z,y))
}
#recode itemtypes from characters to magic numbers
recitt <- function(itemtypes) {
itemtypes[itemtypes=="gaussian"|itemtypes=="normal"|itemtypes=="norm"]=1
itemtypes=as.numeric(itemtypes)
itemtypes
}
# function pp generates vector with parameter names for each itemtype
pp<-function(x){
if(x > 1) return(paste("p",1:x)) # multinomial
if(x== 1) return(c("mean", "stddev")) # normal
}
# returns the number of parameters for each density
np <- function(x) {
if(x > 1) {
return(x)
} else {
if(x==1) {
return(2)
} else {
stop("Unknown itemtype code in function np.")
}
}
}
# convert conpat vector to rows of constraint matrix
pa2conr <- function(x) {
fix=as.logical(x)
x=replace(x,list=which(x==1),0)
un=setdiff(unique(x),0)
y=matrix(0,0,length(x))
for(i in un) {
z=which(x==i)
for(j in 2:length(z)) {
k=rep(0,length(x))
k[z[1]]=1
k[z[j]]=-1
y=rbind(y,k)
}
}
pa = list(fix=fix,conr=y)
return(pa)
}
##################################
# #
# PARAMETER INDEX GENERATION #
# #
##################################
paridx <- function(nstates,itemtypes,mat,idx1=0,idx2=0,it=0,comp=1,group=1) {
# all denotes all parameters of a component model, ie not for a group, and hence
# it does not include the mixing proportions
MATS=c("mix","trans","obser","init","all")
mat = pmatch(mat,MATS)
nitems=length(itemtypes)
nrcomp=length(nstates)
nstates=as.vector(nstates)
## nr of pars per item, the sum of this is the nr of pars per state
# for the response parameters
lobs=sapply(itemtypes,FUN=np)
## nr of pars per component
nrcp = rep(0,nrcomp)
for(i in 1:nrcomp) nrcp[i]=nstates[i]*(nstates[i]+1)+sum(lobs)*nstates[i]
ngpars=nrcomp+sum(nrcp)
mix=1:nrcomp
# mix prop pars
if(mat==1) if(idx1==0) idx=mix else idx=mix[idx1]
# trans pars
if(mat==2) {
tridx <- matrix(1:(nstates[comp]*nstates[comp]),nstates[comp],nstates[comp],byrow=TRUE)
if(idx1==0) idx=c(t(tridx))
else { if(idx2==0) {
idx=tridx[idx1,]
}
else idx=tridx[idx1,idx2]
}
}
# obser pars
if(mat==3) {
en=nstates[comp]*nstates[comp]
obidx=matrix((en+1):(nstates[comp]*sum(lobs)+en),nstates[comp],sum(lobs),byrow=TRUE)
if(it==0) {
if(idx1==0) { idx=c(t(obidx)) }
else {
if(idx2==0) {
idx=obidx[idx1,]
}
else idx=obidx[idx1,idx2]
}
} else { # it>0
if(nitems>1 && it>1) npprev=sum(sapply(itemtypes[1:(it-1)],FUN=np))
else npprev=0
if(itemtypes[it]==1) np=2
else np=itemtypes[it]
if(idx1==0) { #return indices of item parameters for all states
idx=c(t(obidx[,(npprev+1):(npprev+np)]))
} else {
if(idx2==0) idx=obidx[idx1,(npprev+1):(npprev+np)]
else idx=obidx[idx1,npprev+idx2]
}
}
}
# init pars
if(mat==4) {
inidx=(nrcp[comp]-nstates[comp]+1):(nrcp[comp])
if(idx1==0) { idx=inidx }
else idx=inidx[idx1]
}
if(mat==5) idx=1:nrcp[comp]
if(comp>1) idx=idx+sum(nrcp[1:(comp-1)])
if(mat!=1) idx=idx+nrcomp
if(group>1) idx = idx + (group-1)*ngpars
return(idx)
}
# function fresp generates random deviates for itemtype x, given 'pars', ie the item distribution parameters
fresp <- function(x,pars) {
if(x > 1) resp = which(rmultinom(1,size=1,prob=pars)==1) # multinomial
if(x==1) resp = rnorm(1,pars[1],pars[2]) # normal (mu - mean, sigma - stddev)
return(resp)
}
# convert a list of matrices x to a block diagonal matrix
# shamelessly stolen from package assist
bdiag <- function (x) {
if (is.matrix(x))
return(x)
if (!is.list(x))
stop("dismatched input")
n <- length(x)
x <- lapply(x, function(y) if (length(y))
as.matrix(y)
else stop("Zero-length component in x"))
d <- array(unlist(lapply(x, dim)), c(2, n))
rr <- d[1, ]
cc <- d[2, ]
rsum <- sum(rr)
csum <- sum(cc)
out <- array(0, c(rsum, csum))
ind <- array(0, c(4, n))
rcum <- cumsum(rr)
ccum <- cumsum(cc)
ind[1, -1] <- rcum[-n]
ind[2, ] <- rcum
ind[3, -1] <- ccum[-n]
ind[4, ] <- ccum
imat <- array(1:(rsum * csum), c(rsum, csum))
iuse <- apply(ind, 2, function(y, imat) {
imat[(y[1] + 1):y[2], (y[3] + 1):y[4]]
}, imat = imat)
iuse <- as.vector(unlist(iuse))
out[iuse] <- unlist(x)
out
}
#
# STARTING VALUE GENERATION (AUXILIARY) FUNCTIONS
#
kmstart <- function(dat,dmm) {
nst=dmm$nstates
kmc=try(kmeans(na.omit(dat[,1:(nitems(dat)-ncov(dat))]),nst)$cluster,silent=TRUE)
if(is(kmc,"try-error")) {
while(is(kmc,"try-error")) {
kmc=try(kmeans(dat[,1:(nitems(dat)-ncov(dat))],nst)$cluster,silent=TRUE)
}
}
st=cl2st(kmc,dat,dmm)
st
}
poststart <- function(dat,dmm) {
post=posterior(dat,dmm)
st=cl2st(post$states[[1]][,2],dat,dmm)
st
}
cl2st <- function(cluster,dat,dmm) {
nst=dmm$nstates
stin=numeric(nst)
if(is(dmm, "lcm") || dmm$nst==1) {
sttr=diag(nst)
for(i in 1:nst) {
stin[i]=length(cluster[which(cluster==i)])
}
} else {
sttr=matrix(0,nst,nst)
for(i in 1:ind(dat)) {
bg=ifelse(i>1,sum(ntimes(dat)[1:(i-1)])+1,1)
en=bg+ntimes(dat)[i]-1
# add diagonal elements to make sure that the transition table is complete, which
# is not neccessarily so for each case, ie not all transitions necc occur
cl1=c(as.numeric(gl(nst,2)),cluster[bg:en])
cl2=cl1[-1]
cl1=cl1[-length(cl1)]
sttr=sttr+matrix(table(cl1,cl2),nst,nst)
stin[cluster[bg]]=stin[cluster[bg]]+1
# correct for the added transition from nst to the start state
sttr[nst,cluster[bg]]=sttr[nst,cluster[bg]]-1
}
# correct for the added transitions between i and i+1
for(i in 1:nst) {
sttr[i,i]=sttr[i,i]-ind(dat)
if(i<nst) sttr[i,i+1]=sttr[i,i+1]-ind(dat)
}
sttr=sttr/rowSums(sttr)
}
if(any(stin==0)) stin=stin+0.05
stin=stin/sum(stin)
stob=cl2stob(cluster,dat,dmm)
st=c(1,t(sttr),stob,stin)
st
}
tr2stin <- function(sttr) {
etr=eigen(t(sttr))
stin=etr$vectors[,which(etr$values==1)]
stin=stin/sum(stin)
stin
}
cl2stob <- function(cluster,dat,dmm) {
stob=numeric(0)
nst=dmm$nstates
for(i in 1:nst) {
for(j in 1:length(dmm$itemtypes)) {
if(dmm$itemtypes[j]==1) {
meanob=mean(dat[which(cluster==i),j])
sdob=sd(dat[which(cluster==i),j])
stob=c(stob,meanob,sdob)
}
if(dmm$itemtypes[j]>1) {
y=sort(unique(na.omit(dat[,j])))
catprobs=numeric(dmm$itemtypes[j])
for(k in 1:length(y)) {
catprobs[k]=length(dat[which(cluster==i & dat[,j]==y[k]),j])
}
catprobs=catprobs/sum(catprobs)
stob=c(stob,catprobs)
}
}
}
stob
}
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.