[Depmix-commits] r620 - in pkg/depmix: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 20 09:30:40 CET 2014
Author: ingmarvisser
Date: 2014-02-20 09:30:36 +0100 (Thu, 20 Feb 2014)
New Revision: 620
Modified:
pkg/depmix/DESCRIPTION
pkg/depmix/R/depmix-internal.R
pkg/depmix/R/depmix.R
pkg/depmix/R/dmm.R
pkg/depmix/R/generate.R
pkg/depmix/R/lca.R
pkg/depmix/R/mgdmm.R
pkg/depmix/R/mixdmm.R
pkg/depmix/man/fitdmm.Rd
pkg/depmix/src/matrix.cc
pkg/depmix/src/mdmm.cc
pkg/depmix/src/mgdmm.cc
Log:
Resolved issues with global variables and (C/C++) memory issues
Modified: pkg/depmix/DESCRIPTION
===================================================================
--- pkg/depmix/DESCRIPTION 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/DESCRIPTION 2014-02-20 08:30:36 UTC (rev 620)
@@ -1,6 +1,6 @@
Package: depmix
-Version: 0.9.10
-Date: 2011-10-24
+Version: 0.9.11
+Date: 2014-02-19
Title: Dependent Mixture Models
Author: Ingmar Visser <i.visser at uva.nl>
Maintainer: Ingmar Visser <i.visser at uva.nl>
Modified: pkg/depmix/R/depmix-internal.R
===================================================================
--- pkg/depmix/R/depmix-internal.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/depmix-internal.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -44,8 +44,8 @@
dcov=list()
datsplit=list()
for(i in 1:xgmod$ng) {
- dcov[[i]]=markovdata(dat=dat[[i]][,(nrit+1):(dim(dat[[i]])[2]),drop=FALSE],itemt=itt,nt=ntimes(dat[[i]]))
- datsplit[[i]]=markovdata(dat=dat[[i]][,1:nrit,drop=FALSE],itemtypes=itemtypes(dat[[i]])[1:nrit],nt=ntimes(dat[[i]]),replicates=replicates(dat[[i]]))
+ 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.")
Modified: pkg/depmix/R/depmix.R
===================================================================
--- pkg/depmix/R/depmix.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/depmix.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -85,7 +85,7 @@
}
}
- initLogl=loglike(dat=dat,dmm=xgmod,tdcov=tdcov,print=0,set=FALSE)$logl
+ initLogl=loglike(dat=dat,dmm=xgmod,tdcov=tdcov,printlevel=0,set=FALSE)$logl
if(is.nan(initLogl) || is.infinite(initLogl)) stop("Initial model inadmissible, better starting values needed?!\n")
if(printlevel>0) cat("Initial loglikelihood: ", initLogl, " \n")
@@ -149,10 +149,10 @@
# after defining the model
if(nrow(A)>0) {
- Bcomplete=matrix(as.logical(A),nr=nrow(A))
+ Bcomplete=matrix(as.logical(A),nrow=nrow(A))
xcomplete=apply(Bcomplete,1,sum)
A[,which(fixed==0)]=0
- B=matrix(as.logical(A),nr=nrow(A))
+ B=matrix(as.logical(A),nrow=nrow(A))
x=apply(B,1,sum)
for(i in 1:nrow(A)) {
if(x[i]!=xcomplete[i]) {
@@ -277,13 +277,13 @@
# define loglike function
logl <- function(pars) {
xgmod$pars[which(fixed==1)]=pars
- -loglike(dat=dat,dmm=xgmod,print=0,set=FALSE,tdcov=tdcov)$logl
+ -loglike(dat=dat,dmm=xgmod,printlevel=0,set=FALSE,tdcov=tdcov)$logl
}
if(der) {
grad <- function(pars) {
xgmod$pars[which(fixed==1)]=pars
- gr <- -loglike(dat=dat,dmm=xgmod,print=0,set=FALSE,tdcov=tdcov,grad=TRUE)$gr[which(fixed==1)]
+ gr <- -loglike(dat=dat,dmm=xgmod,printlevel=0,set=FALSE,tdcov=tdcov,grad=TRUE)$gr[which(fixed==1)]
return(gr)
}
attr(logl, "gr") <- grad
@@ -297,7 +297,7 @@
lin.upper=bulin,
lin.lower=bllin,
nlin = list(),
- control=donlp2.control(),
+ control=donlp2Control(),
env=.GlobalEnv, name="Rdonlp2")
)
@@ -331,15 +331,15 @@
A=A[,which(constr==1),drop=FALSE]
# remove rows correponding with inequality constraints
idx=which(bllin<bulin)
- ci=matrix(bulin,nc=1)
+ ci=matrix(bulin,ncol=1)
ineq=FALSE
if(length(idx)>0) {
A=A[-idx,,drop=FALSE]
ci=ci[-idx]
Aineq=xgmod$A[idx,which(constr==1),drop=FALSE]
- ciIneqLow=matrix(xgmod$bllin[idx],nc=1)
- ciIneqUp=matrix(xgmod$bulin[idx],nc=1)
+ ciIneqLow=matrix(xgmod$bllin[idx],ncol=1)
+ ciIneqUp=matrix(xgmod$bulin[idx],ncol=1)
ineq=TRUE
}
@@ -375,38 +375,42 @@
violup=violationPenalty*sum(as.numeric(xgpars>xgmod$bu)*sqrt((xgpars-xgmod$bu)^2))
viollow=violationPenalty*sum(as.numeric(xgpars<xgmod$bl)*sqrt((xgpars-xgmod$bl)^2))
violation=violcon+violup+viollow+sum(violineq)
- violation
+ return(violation)
}
- assign("fcalls",0,envir=.GlobalEnv)
- assign("totalviolation",0,envir=.GlobalEnv)
+# assign("fcalls",0,envir=.GlobalEnv)
+# assign("totalviolation",0,envir=.GlobalEnv)
+
+ totalviolation = 0
# optim optimization
if(method=="optim") {
f <- function(pars) {
if(nnewpars>0) xgmod$pars[which(constr==1)]=ginvA%*%ci + b%*%pars[1:nnewpars]
if(nucpars>0) xgmod$pars[which(uncon==1)]=pars[(nnewpars+1):(nnewpars+nucpars)]
- ll=loglike(dat=dat,dmm=xgmod,print=0,set=FALSE,tdcov=tdcov)$logl
+ ll=loglike(dat=dat,dmm=xgmod,printlevel=0,set=FALSE,tdcov=tdcov)$logl
xgpars=xgmod$pars
violation=viol(xgpars)
- assign("totalviolation",violation,envir=.GlobalEnv)
- assign("fcalls",fcalls+1,envir=.GlobalEnv)
+# assign("totalviolation",violation,envir=.GlobalEnv)
+# totalviolation=violation
+ totalviolation=viol(xgpars)
+# assign("fcalls",fcalls+1,envir=.GlobalEnv)
if(is.nan(ll) || is.infinite(ll)) ll=initLogl*1.01
- ll=ll-violation
+ ll=ll-totalviolation
pr=0
- if(printlevel>9 && (fcalls%%10)==0) {cat("call " , fcalls, " logl ", ll, "\n"); pr=1}
- if(printlevel>14 && pr==0) cat("call " , fcalls, " logl ", ll, "\n")
+# if(printlevel>9 && (fcalls%%10)==0) {cat("call " , fcalls, " logl ", ll, "\n"); pr=1}
+# if(printlevel>14 && pr==0) cat("call " , fcalls, " logl ", ll, "\n")
ll
}
gr <- function(pars) {
if(nnewpars>0) xgmod$pars[which(constr==1)]=ginvA%*%ci + b%*%pars[1:nnewpars]
if(nucpars>0) xgmod$pars[which(uncon==1)]=pars[(nnewpars+1):(nnewpars+nucpars)]
- gr=loglike(dat=dat,dmm=xgmod,grad=TRUE,print=0,set=FALSE)$gr
+ gr=loglike(dat=dat,dmm=xgmod,grad=TRUE,printlevel=0,set=FALSE)$gr
if(nnewpars>0) grad=c(gr[which(constr==1)]%*%b,gr[which(uncon==1)])
else grad=gr[which(uncon==1)]
- if(printlevel>0 && (fcalls%%10)==0) cat("call " , fcalls, " grad ", grad, "\n")
- if(printlevel>9) cat("call " , fcalls, " grad ", grad, "\n")
+# if(printlevel>0 && (fcalls%%10)==0) cat("call " , fcalls, " grad ", grad, "\n")
+# if(printlevel>9) cat("call " , fcalls, " grad ", grad, "\n")
grad
}
@@ -444,7 +448,7 @@
# nlm optimization
if(method=="nlm") {
- assign("loglvalue",0,envir=.GlobalEnv)
+# assign("loglvalue",0,envir=.GlobalEnv)
f <- function(pars) {
# expand pars into the model
@@ -452,19 +456,21 @@
if(nucpars>0) xgmod$pars[which(uncon==1)]=pars[(nnewpars+1):(nnewpars+nucpars)]
# get the actual likelihood and gradients if requested/feasible
if(der) {
- llgr=loglike(dat=dat,dmm=xgmod,print=0,grad=TRUE,set=FALSE)
+ llgr=loglike(dat=dat,dmm=xgmod,printlevel=0,grad=TRUE,set=FALSE)
if(nnewpars>0) grad=c(-llgr$gr[which(constr==1)]%*%b,-llgr$gr[which(uncon==1)])
else grad=-llgr$gr[which(uncon==1)]
- } else llgr=loglike(dat=dat,dmm=xgmod,tdcov=tdcov,print=0,grad=FALSE,set=FALSE)
+ } else llgr=loglike(dat=dat,dmm=xgmod,tdcov=tdcov,printlevel=0,grad=FALSE,set=FALSE)
ll=llgr$logl
violation=viol(xgmod$pars)
- assign("totalviolation",violation,envir=.GlobalEnv)
- ll=ll-violation
+ totalviolation=viol(xgmod$pars)
+# totalviolation=violoation
+# assign("totalviolation",violation,envir=.GlobalEnv)
+ ll=ll-totalviolation
pr=0
- assign("fcalls",fcalls+1,envir=.GlobalEnv)
+# assign("fcalls",fcalls+1,envir=.GlobalEnv)
if(is.nan(ll) || is.infinite(ll)) ll=initLogl*1.01
- if(printlevel>9 && (fcalls%%10)==0) {cat("call " , fcalls, " logl ", ll, "\n"); pr=1}
- if(printlevel>14 && pr==0) cat("call " , fcalls, " logl ", ll, "\n")
+# if(printlevel>9 && (fcalls%%10)==0) {cat("call " , fcalls, " logl ", ll, "\n"); pr=1}
+# if(printlevel>14 && pr==0) cat("call " , fcalls, " logl ", ll, "\n")
if(der) attr(ll,'gradient')=grad
-ll
}
@@ -494,7 +500,7 @@
}
names(z)=c("objf","par","grad","inform","iter")
- z$objf=ifelse(loglvalue==0,-z$objf,loglvalue)
+# z$objf=ifelse(loglvalue==0,-z$objf,loglvalue)
} # end nlm
# put fitted pars into xgmod, gather output into z
@@ -687,7 +693,7 @@
post$comp <- list()
for(g in 1:xgmod$ng) {
post$comp[[g]]=matrix(0,length(ntimes(dat[[g]])),xgmod$nrcomp+1)
- post$states[[g]] <- matrix(0,nc=(2+sum(xgmod$nstates)),nr=0)
+ post$states[[g]] <- matrix(0,ncol=(2+sum(xgmod$nstates)),nrow=0)
###the first entry in each series is the number of the component
for(i in 1:(length(ntimes(dat[[g]])))) {
pp=ntimes(dat[[g]])[i]*sum(xgmod$nstates)
@@ -701,7 +707,7 @@
as.integer(i),
PACKAGE="depmix")
post$comp[[g]][i,xgmod$nrcomp+1]=z$postcomp
- pst=matrix(c(rep(z$postcomp,ntimes(dat[[g]])[i]),z$states,matrix(z$postdelta,nc=sum(xgmod$nstates),byrow=TRUE)),nc=(2+sum(xgmod$nstates)))
+ pst=matrix(c(rep(z$postcomp,ntimes(dat[[g]])[i]),z$states,matrix(z$postdelta,ncol=sum(xgmod$nstates),byrow=TRUE)),ncol=(2+sum(xgmod$nstates)))
postprob=pst[ntimes(dat[[g]])[i],3:(2+sum(xgmod$nstates))]
bg=1; en=xgmod$nstates[1];
post$comp[[g]][i,1] = sum(postprob[bg:en])
@@ -724,7 +730,7 @@
cat("Computing standard errors \n")
mod=dmm
eps=1e-8
- lgh=loglike(dat=dat,dmm=mod,grad=TRUE,tdcov=0,print=0)
+ lgh=loglike(dat=dat,dmm=mod,grad=TRUE,tdcov=0,printlevel=0)
ll=lgh$logl
gr=lgh$gr
hsall=matrix(0,mod$npars,mod$npars)
@@ -733,7 +739,7 @@
del=eps*max(1,mod$pars[j])
por=mod$pars[j]
mod$pars[j] = por+del
- alli=loglike(dat=dat,dmm=mod,grad=TRUE,hess=FALSE,set=FALSE,tdcov=0,print=0)
+ alli=loglike(dat=dat,dmm=mod,grad=TRUE,hess=FALSE,set=FALSE,tdcov=0,printlevel=0)
grd=alli$gr
for(i in 1:mod$npars) hsall[i,j] = (grd[i]-gr[i])/del
mod$pars[j]=por
@@ -742,7 +748,7 @@
se=rep(0,mod$npars)
hs=-hsall[which(mod$fixed[1:mod$npars]==1),which(mod$fixed[1:mod$npars]==1)]
if(nrow(mod$A)>0) {
- A=matrix(mod$A[,which(mod$fixed[1:mod$npars]==1)],nr=nrow(mod$A))
+ A=matrix(mod$A[,which(mod$fixed[1:mod$npars]==1)],nrow=nrow(mod$A))
idx=which(mod$bllin<mod$bulin)
if(length(idx)>0) A=A[-idx,,drop=FALSE]
idx=numeric(0)
@@ -799,9 +805,9 @@
PACKAGE="depmix")
}
for(i in 1:samples) {
- gen=generate(dmm=mod, nt=ntimes(dat))
- ll=loglike(dat=gen,dmm=mod,print=0)$logl
- fit=fitdmm(dat=gen,dmm=mod,print=0,ses=0,postst=0,poster=FALSE)
+ gen=generate(dmm=mod, ntimes=ntimes(dat))
+ ll=loglike(dat=gen,dmm=mod,printlevel=0)$logl
+ fit=fitdmm(dat=gen,dmm=mod,printlevel=0,ses=0,postst=0,poster=FALSE)
bootpars[i,]=fit$mod$pars
ll=fit$loglike
gof[i,]=c(ll,fit$aic,fit$bic)
Modified: pkg/depmix/R/dmm.R
===================================================================
--- pkg/depmix/R/dmm.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/dmm.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -51,12 +51,12 @@
# provide parameter start values if none are given, scale start values if provided
mixprop=1.0 # mixing proportion (added for consistency)
if(st) {
- trans <- matrix(stval[2:(nstates*nstates+1)],nr=nstates,byrow=TRUE)
- obser <- matrix(stval[(nstates*nstates+2):(npars-nstates)],nr=nstates,byrow=TRUE)
+ trans <- matrix(stval[2:(nstates*nstates+1)],nrow=nstates,byrow=TRUE)
+ obser <- matrix(stval[(nstates*nstates+2):(npars-nstates)],nrow=nstates,byrow=TRUE)
init=stval[(npars-nstates+1):npars]
} else {
- trans <- matrix(runif(nstates*nstates,0,1),nr=nstates)
- obser <- matrix(0,nr=nstates,nc=sum(lobs))
+ trans <- matrix(runif(nstates*nstates,0,1),nrow=nstates)
+ obser <- matrix(0,nrow=nstates,ncol=sum(lobs))
for(j in 1:nstates) {
z=numeric(0)
for (i in 1:nitems) z=ppar(itemtypes[i],z)
@@ -73,8 +73,8 @@
for(j in 1:nstates) {
for(i in 1:nitems) {
if(itemtypes[i]>1) {
- itpars=pars[paridx(nstates,itemtypes,m="ob",idx1=j,it=i)]
- pars[paridx(nstates,itemtypes,m="ob",idx1=j,it=i)] = itpars/sum(itpars)
+ itpars=pars[paridx(nstates,itemtypes,mat="ob",idx1=j,it=i)]
+ pars[paridx(nstates,itemtypes,mat="ob",idx1=j,it=i)] = itpars/sum(itpars)
}
}
}
@@ -102,7 +102,7 @@
blst=numeric(0); bust=numeric(0)
if(nstates>1) {
sctr=matrix(0,nstates,npars)
- for(i in 1:nstates) sctr[i,paridx(nstates,itemtypes,m="tr",i)]=1
+ for(i in 1:nstates) sctr[i,paridx(nstates,itemtypes,mat="tr",i)]=1
blst <- rep(1,nstates)
bust <- rep(1,nstates)
}
@@ -112,7 +112,7 @@
for(j in 1:nitems) {
if(itemtypes[j]>1) {
sc=rep(0,npars)
- sc[paridx(nstates,itemtypes,m="ob",it=j,idx1=i)]=1
+ sc[paridx(nstates,itemtypes,mat="ob",it=j,idx1=i)]=1
scob=rbind(scob,sc)
blso=c(blso,1); buso=c(buso,1)
}
@@ -133,7 +133,7 @@
else A=rbind(sctr,scob,scin)
# add user supplied linear constraints to A and set their bounds
if(cr) {
- conrows <- matrix(conrows,nc=npars,byrow=TRUE)
+ conrows <- matrix(conrows,ncol=npars,byrow=TRUE)
A=rbind(A,conrows)
bllin <- c(bllin,rep(0,nrow(conrows)))
bulin <- c(bulin,rep(0,nrow(conrows)))
@@ -152,23 +152,23 @@
td=1 # only one covariate supported at the moment
bltd=rep(0,npars)
butd=rep(0,npars)
- if(any(tdfix[paridx(nstates,itemtypes,m="tr")]==1)) {
+ if(any(tdfix[paridx(nstates,itemtypes,mat="tr")]==1)) {
tdtr=TRUE
if(nstates==1) stop("Covariates on the transition matrix do not make sense with a single state model.\n")
- butd[paridx(nstates,itemtypes,m="tr")]=1
- bltd[paridx(nstates,itemtypes,m="tr")]=-1
+ butd[paridx(nstates,itemtypes,mat="tr")]=1
+ bltd[paridx(nstates,itemtypes,mat="tr")]=-1
}
# etc for other betas
- if(any(tdfix[paridx(nstates,itemtypes,m="ob")]==1)) {
+ if(any(tdfix[paridx(nstates,itemtypes,mat="ob")]==1)) {
tdob=TRUE
- butd[paridx(nstates,itemtypes,m="ob")]=bigB
- bltd[paridx(nstates,itemtypes,m="ob")]=-bigB
+ butd[paridx(nstates,itemtypes,mat="ob")]=bigB
+ bltd[paridx(nstates,itemtypes,mat="ob")]=-bigB
}
- if(any(tdfix[paridx(nstates,itemtypes,m="in")]==1)) {
+ if(any(tdfix[paridx(nstates,itemtypes,mat="in")]==1)) {
tdin=TRUE
if(nstates==1) stop("Covariates on the initial probs do not make sense with a single state model.\n")
- butd[paridx(nstates,itemtypes,m="in")]=1
- bltd[paridx(nstates,itemtypes,m="in")]=-1
+ butd[paridx(nstates,itemtypes,mat="in")]=1
+ bltd[paridx(nstates,itemtypes,mat="in")]=-1
}
tdpars=rep(0,npars)
if(is.null(tdst)) tdpars=replace(tdpars,tdfix==1,rnorm(sum(tdfix),0,0.1))
@@ -177,20 +177,20 @@
# covariates for tr parameters
if(tdtr) {
# values: make values conformable with constraints
- tdtrans=matrix(tdpars[paridx(nstates,itemtypes,m="tr")],nr=nstates,byrow=TRUE)
+ tdtrans=matrix(tdpars[paridx(nstates,itemtypes,mat="tr")],nrow=nstates,byrow=TRUE)
if(nstates==2) tdtrans[,2]=-tdtrans[,1]
else tdtrans[,nstates]=-apply(tdtrans[,(1:(nstates-1))],1,sum)
- tdpars[paridx(nstates,itemtypes,m="tr")]=t(tdtrans)
+ tdpars[paridx(nstates,itemtypes,mat="tr")]=t(tdtrans)
# make the constraints
sctdtr=matrix(0,nstates+nstates*nstates,2*npars)
blsctdtr=rep(0,nstates+nstates*nstates)
busctdtr=c(rep(0,nstates),rep(1,nstates*nstates))
# ... between the beta's (these sum to zero per state)
- for(i in 1:nstates) sctdtr[i,paridx(nstates,itemtypes,m="tr",idx1=i)+npars]=1
+ for(i in 1:nstates) sctdtr[i,paridx(nstates,itemtypes,mat="tr",idx1=i)+npars]=1
# ... between transpars and beta's 0<=a_ij+b_ij<=1, assuming the covariate is between 0 and 1
for(i in 1:nstates) {
for(j in 1:nstates) {
- pidx=paridx(nstates,itemtypes,m="tr",idx1=i,idx2=j)
+ pidx=paridx(nstates,itemtypes,mat="tr",idx1=i,idx2=j)
sctdtr[nstates+(i-1)*nstates+j,c(pidx,pidx+npars)]=c(1,1)
}
}
@@ -208,10 +208,10 @@
if(itemtypes[j]>1) {
kk=kk+1
# these beta's have to sum to one
- pp=tdpars[paridx(nstates,itemtypes,m="ob",idx1=i,it=j)]
+ pp=tdpars[paridx(nstates,itemtypes,mat="ob",idx1=i,it=j)]
pp[itemtypes[j]]=-sum(pp[1:(itemtypes[j]-1)])
- tdpars[paridx(nstates,itemtypes,m="ob",idx1=i,it=j)]=pp
- sctdob[kk,paridx(nstates,itemtypes,m="ob",idx1=i,it=j)+npars]=1
+ tdpars[paridx(nstates,itemtypes,mat="ob",idx1=i,it=j)]=pp
+ sctdob[kk,paridx(nstates,itemtypes,mat="ob",idx1=i,it=j)+npars]=1
}
}
}
@@ -222,7 +222,7 @@
if(itemtypes[j]>1) {
for(k in 1:itemtypes[j]) {
kk=kk+1
- pidx=paridx(nstates,itemtypes,m="ob",idx1=i,it=j,idx2=k)
+ pidx=paridx(nstates,itemtypes,mat="ob",idx1=i,it=j,idx2=k)
sctdob[kk,c(pidx,pidx+npars)]=c(1,1)
}
}
@@ -240,10 +240,10 @@
blsctdin=rep(0,1+nstates)
busctdin=c(0,rep(1,nstates))
# beta's sum to zero
- sctdin[,paridx(nstates,itemtypes,m="in")+npars]=1
+ sctdin[,paridx(nstates,itemtypes,mat="in")+npars]=1
# beta + initpar should be between zero and one
for(i in 1:nstates) {
- pidx=paridx(nstates,itemtypes,m="in",idx1=i)
+ pidx=paridx(nstates,itemtypes,mat="in",idx1=i)
sctdin[1+i,c(pidx,pidx+npars)]=c(1,1)
}
}
@@ -281,10 +281,10 @@
# check for zero rows in A which may result from fixing parameters at zero, and remove them
# set zeroes in A when pars are fixed, and change bllin and bulin accordingly
if(nrow(A)>0) {
- Bcomplete=matrix(as.logical(A),nr=nrow(A))
+ Bcomplete=matrix(as.logical(A),nrow=nrow(A))
xcomplete=apply(Bcomplete,1,sum)
A[,which(fixed==0)]=0
- B=matrix(as.logical(A),nr=nrow(A))
+ B=matrix(as.logical(A),nrow=nrow(A))
x=apply(B,1,sum)
for(i in 1:nrow(A)) {
if(x[i]!=xcomplete[i]) {
@@ -330,16 +330,16 @@
ses=ifelse(is.null(se),0,1)
if(!"lcm" %in% class(object)) {
cat(" Parameter values, transition matrix \n\n")
- trsx=matrix(object$pars[paridx(object$nstates,object$itemtypes,m="tr")],object$nstates,byrow=TRUE)
+ trsx=matrix(object$pars[paridx(object$nstates,object$itemtypes,mat="tr")],object$nstates,byrow=TRUE)
rownames(trsx)=object$snames
if(object$tdtr && object$tdfit==1) {
- betr=matrix(object$pars[paridx(object$nstates,object$itemtypes,m="tr")+object$npars],object$nstates,byrow=TRUE)
+ betr=matrix(object$pars[paridx(object$nstates,object$itemtypes,mat="tr")+object$npars],object$nstates,byrow=TRUE)
trsx=matrix(rbind(c(trsx),c(betr)),object$nstates*2,object$nstates)
be=rep(x="be", object$nstates)
rownames(trsx)=c(rbind(object$snames,be))
}
if(ses) {
- setr=matrix(se[paridx(object$nstates,object$itemtypes,m="tr")],object$nstates,byrow=TRUE)
+ setr=matrix(se[paridx(object$nstates,object$itemtypes,mat="tr")],object$nstates,byrow=TRUE)
tr=matrix(0,object$nstates*3,object$nstates)
for(i in 1: object$nstates) {
tr[3*(i-1)+1,]=trsx[i,]
@@ -356,17 +356,17 @@
cat("\n\n")
}
cat(" Parameter values, observation parameters \n\n")
- obsx=matrix(object$pars[paridx(object$nstates,object$itemtypes,m="ob")],object$nstates,byrow=TRUE)
+ obsx=matrix(object$pars[paridx(object$nstates,object$itemtypes,mat="ob")],object$nstates,byrow=TRUE)
rownames(obsx)=object$snames
if(object$tdob && object$tdfit==1) {
- tdobpars=matrix(object$pars[paridx(object$nstates,object$itemtypes,m="ob")+object$npars],object$nstates,byrow=TRUE)
+ tdobpars=matrix(object$pars[paridx(object$nstates,object$itemtypes,mat="ob")+object$npars],object$nstates,byrow=TRUE)
obsx=matrix(rbind(c(obsx),c(tdobpars)),object$nstates*2,byrow=FALSE)
be=rep(x="be", object$nstates)
rownames(obsx)=c(rbind(object$snames,be))
}
if(ses) {
- seob=matrix(se[paridx(object$nstates,object$itemtypes,m="ob")],object$nstates,byrow=TRUE)
- ob=matrix(0,object$nstates*3,nc=ncol(obsx))
+ seob=matrix(se[paridx(object$nstates,object$itemtypes,mat="ob")],object$nstates,byrow=TRUE)
+ ob=matrix(0,object$nstates*3,ncol=ncol(obsx))
for(i in 1: object$nstates) {
ob[3*(i-1)+1,]=obsx[i,]
ob[3*(i-1)+2,]=seob[i,]
@@ -387,18 +387,18 @@
if(!"lcm" %in% class(object)) cat(" Parameter values, initial state probabilies \n\n")
else cat(" Parameter values, unconditional (class) probabilities \n\n")
if(ses) {
- inix=matrix(c(object$pars[paridx(object$nstates,object$itemtypes,m="in")],rep(0,2*object$nstates)),3,object$nstates,byrow=TRUE)
- sein=se[paridx(object$nstates,object$itemtypes,m="in")]
+ inix=matrix(c(object$pars[paridx(object$nstates,object$itemtypes,mat="in")],rep(0,2*object$nstates)),3,object$nstates,byrow=TRUE)
+ sein=se[paridx(object$nstates,object$itemtypes,mat="in")]
for(i in 1: object$nstates) {
inix[2,i]=sein[i]
inix[3,i]=ifelse(sein[i]==0,NA,inix[1,i]/sein[i])
}
rownames(inix)=c("val","se","t")
} else {
- inix=matrix(object$pars[paridx(object$nstates,object$itemtypes,m="in")],1,object$nstates,byrow=TRUE)
+ inix=matrix(object$pars[paridx(object$nstates,object$itemtypes,mat="in")],1,object$nstates,byrow=TRUE)
rownames(inix)="val"
if(object$tdin && object$tdfit==1) {
- inix=rbind(inix,object$pars[paridx(object$nstates,object$itemtypes,m="in")+object$npars])
+ inix=rbind(inix,object$pars[paridx(object$nstates,object$itemtypes,mat="in")+object$npars])
rownames(inix)=c("val","be")
}
}
Modified: pkg/depmix/R/generate.R
===================================================================
--- pkg/depmix/R/generate.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/generate.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -11,17 +11,17 @@
trans=list(); obser=list(); init=list()
if(class(dmm)[1]=="dmm") dmm=mixdmm(dmm=list(dmm))
for(i in 1:dmm$nrcomp) {
- idx=paridx(dmm$nstates,dmm$itemtypes,m="tr",comp=i)
+ idx=paridx(dmm$nstates,dmm$itemtypes,mat="tr",comp=i)
trans[[i]]=matrix(dmm$pars[idx],dmm$nstates[i],byrow=TRUE)
- idx=paridx(dmm$nstates,dmm$itemtypes,m="ob",comp=i)
+ idx=paridx(dmm$nstates,dmm$itemtypes,mat="ob",comp=i)
obser[[i]]=matrix(dmm$pars[idx],dmm$nstates[i],byrow=TRUE)
- idx=paridx(dmm$nstates,dmm$itemtypes,m="in",comp=i)
+ idx=paridx(dmm$nstates,dmm$itemtypes,mat="in",comp=i)
init[[i]]=matrix(dmm$pars[idx],dmm$nstates[i])
}
ntcount=0
# create return value vectors
nitems=length(dmm$itemtypes)
- obs = matrix(0,nc=nitems,nr=sum(ntimes))
+ obs = matrix(0,ncol=nitems,nrow=sum(ntimes))
instates=integer(length(ntimes))
incomp=integer(length(ntimes))
for(j in 1:length(ntimes)) {
@@ -48,7 +48,7 @@
ntcount = ntcount + ntimes[j]
}
# repeat ntimes
- gen <- markovdata(dat=obs,itemt=dmm$itemtnames,ntimes=ntimes)
+ gen <- markovdata(dat=obs,itemtypes=dmm$itemtnames,ntimes=ntimes)
attr(gen,"instates")=instates
if(dmm$nrcomp>1) attr(gen,"incomp")=incomp
gen
Modified: pkg/depmix/R/lca.R
===================================================================
--- pkg/depmix/R/lca.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/lca.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -32,13 +32,13 @@
stval=c(1,diag(nclasses),runif((lcmnp-1)))
}
if(!is.null(conrows)) {
- conrows=matrix(conrows,nc=lcmnp,byrow=TRUE)
- trzeroes=matrix(0,nr=nrow(conrows),nc=nclasses*nclasses)
+ conrows=matrix(conrows,ncol=lcmnp,byrow=TRUE)
+ trzeroes=matrix(0,nrow=nrow(conrows),ncol=nclasses*nclasses)
mpzeroes=matrix(0,nrow(conrows),1)
conrows=t(cbind(mpzeroes,trzeroes,conrows[,2:(ncol(conrows))]))
}
if(!is.null(linmat)) {
- trzeroes=matrix(0,nr=nrow(linmat),nc=nclasses*nclasses)
+ trzeroes=matrix(0,nrow=nrow(linmat),ncol=nclasses*nclasses)
linmat=cbind(trzeroes,linmat)
}
if(!is.null(conpat)) {
Modified: pkg/depmix/R/mgdmm.R
===================================================================
--- pkg/depmix/R/mgdmm.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/mgdmm.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -60,11 +60,11 @@
# count the number of transition pars and add this many rows *(ng-1) to A
trcount=0
for(i in 1:nrcomp) {
- trcount = length(paridx(nstates,itemtypes,comp=i,m="tr"))
+ trcount = length(paridx(nstates,itemtypes,comp=i,mat="tr"))
trinv = matrix(0,trcount*(ng-1),npars)
j=0
for(g in 1:(ng-1)) {
- idx = paridx(nstates,itemtypes,comp=i,m="tr")
+ idx = paridx(nstates,itemtypes,comp=i,mat="tr")
for (id in idx) {
j = j+1
trinv[j,id] = 1
@@ -83,11 +83,11 @@
# count the number of observation pars and add this many rows *(ng-1) to A
obcount=0
for(i in 1:nrcomp) {
- obcount = length(paridx(nstates,itemtypes,comp=i,m="ob"))
+ obcount = length(paridx(nstates,itemtypes,comp=i,mat="ob"))
obinv = matrix(0,obcount*(ng-1),npars)
j=0
for(g in 1:(ng-1)) {
- idx = paridx(nstates,itemtypes,comp=i,m="ob")
+ idx = paridx(nstates,itemtypes,comp=i,mat="ob")
for(id in idx) {
j = j+1
obinv[j,id] = 1
@@ -106,11 +106,11 @@
# count the number of initial state pars and add this many rows *(ng-1) to A
incount=0
for(i in 1:nrcomp) {
- incount = length(paridx(nstates,itemtypes,comp=i,m="in"))
+ incount = length(paridx(nstates,itemtypes,comp=i,mat="in"))
ininv = matrix(0,incount*(ng-1),npars)
j=0
for(g in 1:(ng-1)) {
- idx = paridx(nstates,itemtypes,comp=i,m="in")
+ idx = paridx(nstates,itemtypes,comp=i,mat="in")
for (id in idx) {
j = j+1
ininv[j,id] = 1
@@ -180,7 +180,7 @@
idx=numeric(0)
if(nrow(A)>0) {
for(i in 1:nrow(A)) {if(sum(as.logical(A[i,]))!=0) idx = c(idx,i) }
- A=matrix(A[idx,],nc=nparstotal)
+ A=matrix(A[idx,],ncol=nparstotal)
}
# ... also remove corresponding entries in bllin and bulin
bllin=bllin[idx]
Modified: pkg/depmix/R/mixdmm.R
===================================================================
--- pkg/depmix/R/mixdmm.R 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/R/mixdmm.R 2014-02-20 08:30:36 UTC (rev 620)
@@ -48,13 +48,13 @@
fixed=c(fixed,mod[[i]]$fixed[2:mod[[i]]$npars])
}
# sum constraint for mixprop
- sc=matrix(rep(1,nrcomp),nr=1)
+ sc=matrix(rep(1,nrcomp),nrow=1)
bllin=1
bulin=1
# user supplied general linear constraints on mixprop
amp=NULL
if(!(is.null(conrows))) {
- conrows=matrix(conrows,nc=npars,byrow=TRUE)
+ conrows=matrix(conrows,ncol=npars,byrow=TRUE)
amp=conrows
bllin=c(bllin,rep(0,nrow(amp)))
bulin=c(bulin,rep(0,nrow(amp)))
@@ -82,7 +82,7 @@
}
}
# make combined constraint matrix
- A=matrix(0,nc=npars,nr=nr)
+ A=matrix(0,ncol=npars,nrow=nr)
A[1:nrow(amp),1:nrcomp] <- amp
B=bdiag(linmat)
A[2:nrow(A),(nrcomp+1):npars]=B
@@ -119,7 +119,7 @@
linmattd[[i]]=matrix(0,1,mod[[i]]$npars-1)
}
}
- zeroes=matrix(0,nc=npars,nr=nr)
+ zeroes=matrix(0,ncol=npars,nrow=nr)
A=cbind(A,zeroes)
C=bdiag(linmattd)
A[2:nrow(A),(npars+nrcomp+1):(2*npars)]=C
Modified: pkg/depmix/man/fitdmm.Rd
===================================================================
--- pkg/depmix/man/fitdmm.Rd 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/man/fitdmm.Rd 2014-02-20 08:30:36 UTC (rev 620)
@@ -259,7 +259,8 @@
stcov=rep(0,15)
stcov[2:5]=c(-0.4,0.4,0.15,-0.15)
-mod<-dmm(nstates=2,itemt=c("n",2),stval=stv,conpat=conpat,tdfix=tdfix,tdst=stcov,modname="twoboth+cov")
+mod<-dmm(nstates=2,itemt=c("n",2),stval=stv,conpat=conpat,tdfix=tdfix,tdst=stcov,
+modname="twoboth+cov")
fit3 <- fitdmm(dat=speed,dmm=mod,tdcov=1,der=0,ses=0,vfa=80)
summary(fit3)
Modified: pkg/depmix/src/matrix.cc
===================================================================
--- pkg/depmix/src/matrix.cc 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/src/matrix.cc 2014-02-20 08:30:36 UTC (rev 620)
@@ -73,20 +73,26 @@
//vector access
double& matrix::operator() (const int x) {
-#ifdef MATRIXBOUNDS
- if(row==1) {
- if(x<1 || x>col) error("[Matrix] Error: rowvector out of range.\n");
- }
- else {
- if(col==1) {
- if (x<1 || x>row) error("[Matrix] Error: colvector out of range.\n");
- }
- else error("[Matrix] Error: matrix adressed as vector.\n");
- }
-#endif
+ #ifdef MATRIXBOUNDS
+ if(row==1) {
+ if(x<1 || x>col) {
+ error("[Matrix] Error: rowvector out of range.\n");
+ }
+ }
+ else {
+ if(col==1) {
+ if (x<1 || x>row) {
+ error("[Matrix] Error: colvector out of range.\n");
+ }
+ } else {
+ error("[Matrix] Error: matrix adressed as vector.\n");
+ }
+ }
+ #endif
if (row==1) return(value[0][x-1]);
else {
if (col==1) return(value[x-1][0]);
+ else error("[Matrix] Error: matrix adressed as vector.\n");
}
}
Modified: pkg/depmix/src/mdmm.cc
===================================================================
--- pkg/depmix/src/mdmm.cc 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/src/mdmm.cc 2014-02-20 08:30:36 UTC (rev 620)
@@ -268,6 +268,7 @@
if(npit!=it) return(0.0);
else {
int npitpar=itempar(np);
+
//multinomial item
if(getMode(it)>1) {
if(npitpar!=((int) datit)) return(0.0);
@@ -283,6 +284,7 @@
}
}
}
+ return(0.0); // to make sure we do not get 'control may reach end of non-void function' warnings
}
}
@@ -312,6 +314,7 @@
}
}
}
+ return(0.0); // to make sure we do not get 'control may reach end of non-void function' warnings
}
}
Modified: pkg/depmix/src/mgdmm.cc
===================================================================
--- pkg/depmix/src/mgdmm.cc 2014-02-05 21:29:21 UTC (rev 619)
+++ pkg/depmix/src/mgdmm.cc 2014-02-20 08:30:36 UTC (rev 620)
@@ -79,7 +79,6 @@
}
mgdmm::~mgdmm(void) {
- delete [] nstates;
delete [] itemtypes;
delete [] mods;
delete [] ncpars;
@@ -97,10 +96,10 @@
delete [] omegat;
delete [] omegafinal;
}
+ delete [] nstates;
}
void mgdmm::reset(const int ngroups, const int nrc, int *nst, const int nit, int *itt, int xm) {
- delete [] nstates;
delete [] itemtypes;
delete [] mods;
delete [] ncpars;
@@ -118,6 +117,7 @@
delete [] omegat;
delete [] omegafinal;
}
+ delete [] nstates;
initialize(ngroups,nrc,nst,nit,itt,xm);
}
More information about the depmix-commits
mailing list