[Depmix-commits] r394 - in trunk: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 8 19:47:27 CET 2010
Author: ingmarvisser
Date: 2010-03-08 19:47:27 +0100 (Mon, 08 Mar 2010)
New Revision: 394
Modified:
trunk/NAMESPACE
trunk/R/depmixfit.R
trunk/R/transInit.R
trunk/depmixNew-test6.R
trunk/man/depmix.fit.Rd
Log:
Minor changes, some resolved conflicts
Modified: trunk/NAMESPACE
===================================================================
--- trunk/NAMESPACE 2010-03-08 16:09:22 UTC (rev 393)
+++ trunk/NAMESPACE 2010-03-08 18:47:27 UTC (rev 394)
@@ -44,7 +44,8 @@
nstates,
depmix,
mix,
- posterior,
+ posterior,
+ getModel,
GLMresponse,
MVNresponse,
transInit,
Modified: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2010-03-08 16:09:22 UTC (rev 393)
+++ trunk/R/depmixfit.R 2010-03-08 18:47:27 UTC (rev 394)
@@ -26,7 +26,6 @@
}
}
-
# check feasibility of starting values
if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is NaN; please provide better starting values. ")
@@ -36,7 +35,10 @@
if(method=="donlp"||method=="rsolnp") {
- if(method=="donlp"&!require(Rdonlp2)) method="rsolnp"
+ if(method=="donlp"&!require(Rdonlp2)) {
+ warning("Rdonlp2 not available, method changed to rsolnp")
+ method="rsolnp"
+ }
if(method=="rsolnp"&!(require(Rsolnp))) stop("Optimization requires either 'Rdonlp2' or 'Rsolnp'")
Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R 2010-03-08 16:09:22 UTC (rev 393)
+++ trunk/R/transInit.R 2010-03-08 18:47:27 UTC (rev 394)
@@ -1,212 +1,212 @@
-#
-# for the transition models and the prior (y is missing, ie there is no
-# response, and nstates must be provided as the number of categories
-# neccessary in the mulinomial model)
-#
-
-setClass("transInit",contains="GLMresponse")
-
-setMethod("transInit",
- signature(formula="formula"),
- function(formula,nstates,data=NULL,family=multinomial(),pstart=NULL,fixed=NULL,prob=TRUE, ...) {
- call <- match.call()
- if(formula==formula(~1) & is.null(data)) {
- x <- matrix(1,ncol=1)
- } else {
- mf <- match.call(expand.dots = FALSE)
- m <- match(c("formula", "data"), names(mf), 0)
- mf <- mf[c(1, m)]
- mf$drop.unused.levels <- TRUE
- mf[[1]] <- as.name("model.frame")
- mf <- eval(mf, parent.frame())
- x <- model.matrix(attr(mf, "terms"),mf)
- }
- y <- matrix(1,ncol=1) # y is not needed in the transition and init models
- parameters <- list()
- constr <- NULL
- if(is.null(nstates)) stop("'nstates' must be provided in call to transInit model")
- if(family$family=="multinomial") {
- if(family$link=="identity") {
- parameters$coefficients <- t(apply(matrix(1,ncol=nstates,nrow=ncol(x)),1,function(x) x/sum(x)))
- if(is.null(fixed)) {
- fixed <- matrix(0,nrow=nrow(parameters$coefficients),ncol=ncol(parameters$coefficients))
- fixed <- rep(0,nstates) # this needs to be fixed at some point using contraints
- fixed <- c(as.logical(fixed))
- }
- constr <- list(
- lin = matrix(1,nr=1,nc=nstates),
- linup = 1,
- linlow = 1,
- parup = rep(1,nstates),
- parlow = rep(0,nstates)
- )
- } else {
- parameters$coefficients <- matrix(0,ncol=nstates,nrow=ncol(x))
- if(is.null(fixed)) {
- fixed <- parameters$coefficients
- fixed[,family$base] <- 1
- fixed <- c(as.logical(t(fixed)))
- }
- }
- }
- npar <- length(unlist(parameters))
- if(is.null(fixed)) fixed <- rep(0,npar)
- if(!is.null(pstart)) {
- if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
- if(family$family=="multinomial") {
- if(family$link=="identity") {
- parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
- parameters$coefficients[1,] <- parameters$coefficients[1,]/sum(parameters$coefficients[1,])
- pstart <- matrix(pstart,ncol(x),byrow=TRUE)
- if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),] # this cannot occur ...
- } else {
- if(prob) {
- parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
- } else {
- parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
- }
- pstart <- matrix(pstart,ncol(x),byrow=TRUE)
- if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
- }
- } else {
- if(family$link=="identity") parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)])
- else parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)],base=family$base)
- }
- }
- # FIX this: do we need a switch here?
- mod <- switch(family$family,
- multinomial = new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr)
- )
- mod
- }
-)
-
-setMethod("logDens","transInit",
- function(object) {
- log(predict(object))
- }
-)
-
-setMethod("dens","transInit",
- function(object,log=FALSE) {
- if(log) log(predict(object))
- else predict(object)
- }
-)
-
-setMethod("predict","transInit",
- function(object,newx) {
- if(missing(newx)) {
- if(object at family$link=="identity") object at family$linkinv(object at x%*%object at parameters$coefficients)
- else object at family$linkinv(object at x%*%object at parameters$coefficients,base=object at family$base)
- } else {
- if(!(is.matrix(newx))) stop("'newx' must be matrix in predict(transInit)")
- if(!(ncol(newx)==nrow(object at parameters$coefficients))) stop("Incorrect dimension of 'newx' in predict(transInit)")
- if(object at family$link=="identity") object at family$linkinv(newx%*%object at parameters$coefficients)
- else object at family$linkinv(newx%*%object at parameters$coefficients,base=object at family$base)
- }
- }
-)
-
-setMethod("fit","transInit",
- function(object,w,ntimes) {
- pars <- object at parameters
- if(missing(w)) w <- NULL
- pars <- object at parameters
- base <- object at family$base # delete me
- #y <- object at y[,-base]
- y <- object at y
- x <- object at x
- if(is.matrix(y)) na <- unlist(apply(y,2,function(x) which(is.na(x)))) else na <- which(is.na(y))
- if(is.matrix(x)) na <- c(na,unlist(apply(x,2,function(x) which(is.na(x))))) else na <- c(na,which(is.na(x)))
- if(!is.null(w)) na <- c(na,which(is.na(w)))
- y <- as.matrix(y)
- x <- as.matrix(x)
- na <- unique(na)
- if(length(na)>0) {
- x <- x[-na,]
- y <- y[-na,]
- #y <- round(y) # delete me
- if(!is.null(w)) w <- w[-na]
- }
- switch(object at family$link,
- mlogit = {
- mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
- mask[,base] <- 0 # fix base category coefficients to 0
- mask <- rbind(0,mask) # fix "bias" nodes to 0
- Wts <- mask
- Wts[-1,] <- pars$coefficients # set starting weights
- Wts[Wts == Inf] <- .Machine$double.max.exp # Fix this!!!!
- Wts[Wts == -Inf] <- .Machine$double.min.exp # Fix this!!!!!
- if(!is.null(w)) {
- if(NCOL(y) < 3) {
- fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- } else {
- fit <- nnet.default(x,y,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- }
- } else {
- if(NCOL(y) < 3) {
- fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- } else {
- fit <- nnet.default(x,y,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- }
- }
- pars$coefficients <- t(matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1)[-1,])
- object <- setpars(object,unlist(pars))
- },
- identity = {
- # object at y = fbo$xi[,,i]/fbo$gamma[,i]
- # should become (see em):
- #for(k in 1:ns) {
- # trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
- # }
- if(!is.null(w)) {
- sw <- sum(w)
- pars <- colSums(w*object at y)/sum(w)
- } else {
- pars <- colMeans(object at y)
- }
- object <- setpars(object,pars)
- },
- stop("link function not implemented")
- )
- object
- }
-)
-
-setMethod("simulate",signature(object="transInit"),
- function(object,nsim=1,seed=NULL,times,is.prior=FALSE,...) {
- if(!is.null(seed)) set.seed(seed)
- if(is.prior) {
- pr <- dens(object)
- sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nrow(pr)))
- states <- t(apply(sims,c(2,3), function(x) which(x==1)))
- return(states)
- } else {
- if(missing(times)) {
- # this is likely to be a stationary model...???
- pr <- predict(object)
- } else {
- pr <- predict(object)[times,]
- if(length(times)==1) pr <- matrix(pr,ncol=length(pr))
- }
- nt <- nrow(pr)
- sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
- states <- t(apply(sims,c(2,3), function(x) which(x==1)))
- # states <- apply(apply(pr,2,rmultinom rmultinom(nt*nsim,size=1,prob=pr),2,function(x) which(x==1))
- return(states)
- }
- }
-)
-
-setMethod("getdf","transInit",
- function(object) {
- npar <- sum(!object at fixed)
- if(object at family$link == "identity") {
- npar <- npar-1
- if(npar<0) npar <- 0
- }
- npar
- }
-)
+#
+# for the transition models and the prior (y is missing, ie there is no
+# response, and nstates must be provided as the number of categories
+# neccessary in the mulinomial model)
+#
+
+setClass("transInit",contains="GLMresponse")
+
+setMethod("transInit",
+ signature(formula="formula"),
+ function(formula,nstates,data=NULL,family=multinomial(),pstart=NULL,fixed=NULL,prob=TRUE, ...) {
+ call <- match.call()
+ if(formula==formula(~1) & is.null(data)) {
+ x <- matrix(1,ncol=1)
+ } else {
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data"), names(mf), 0)
+ mf <- mf[c(1, m)]
+ mf$drop.unused.levels <- TRUE
+ mf[[1]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+ x <- model.matrix(attr(mf, "terms"),mf)
+ }
+ y <- matrix(1,ncol=1) # y is not needed in the transition and init models
+ parameters <- list()
+ constr <- NULL
+ if(is.null(nstates)) stop("'nstates' must be provided in call to transInit model")
+ if(family$family=="multinomial") {
+ if(family$link=="identity") {
+ parameters$coefficients <- t(apply(matrix(1,ncol=nstates,nrow=ncol(x)),1,function(x) x/sum(x)))
+ if(is.null(fixed)) {
+ fixed <- matrix(0,nrow=nrow(parameters$coefficients),ncol=ncol(parameters$coefficients))
+ fixed <- rep(0,nstates) # this needs to be fixed at some point using contraints
+ fixed <- c(as.logical(fixed))
+ }
+ constr <- list(
+ lin = matrix(1,nr=1,nc=nstates),
+ linup = 1,
+ linlow = 1,
+ parup = rep(1,nstates),
+ parlow = rep(0,nstates)
+ )
+ } else {
+ parameters$coefficients <- matrix(0,ncol=nstates,nrow=ncol(x))
+ if(is.null(fixed)) {
+ fixed <- parameters$coefficients
+ fixed[,family$base] <- 1
+ fixed <- c(as.logical(t(fixed)))
+ }
+ }
+ }
+ npar <- length(unlist(parameters))
+ if(is.null(fixed)) fixed <- rep(0,npar)
+ if(!is.null(pstart)) {
+ if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
+ if(family$family=="multinomial") {
+ if(family$link=="identity") {
+ parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+ parameters$coefficients[1,] <- parameters$coefficients[1,]/sum(parameters$coefficients[1,])
+ pstart <- matrix(pstart,ncol(x),byrow=TRUE)
+ if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),] # this cannot occur ...
+ } else {
+ if(prob) {
+ parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
+ } else {
+ parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+ }
+ pstart <- matrix(pstart,ncol(x),byrow=TRUE)
+ if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
+ }
+ } else {
+ if(family$link=="identity") parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)])
+ else parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)],base=family$base)
+ }
+ }
+ # FIX this: do we need a switch here?
+ mod <- switch(family$family,
+ multinomial = new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr)
+ )
+ mod
+ }
+)
+
+setMethod("logDens","transInit",
+ function(object) {
+ log(predict(object))
+ }
+)
+
+setMethod("dens","transInit",
+ function(object,log=FALSE) {
+ if(log) log(predict(object))
+ else predict(object)
+ }
+)
+
+setMethod("predict","transInit",
+ function(object,newx) {
+ if(missing(newx)) {
+ if(object at family$link=="identity") object at family$linkinv(object at x%*%object at parameters$coefficients)
+ else object at family$linkinv(object at x%*%object at parameters$coefficients,base=object at family$base)
+ } else {
+ if(!(is.matrix(newx))) stop("'newx' must be matrix in predict(transInit)")
+ if(!(ncol(newx)==nrow(object at parameters$coefficients))) stop("Incorrect dimension of 'newx' in predict(transInit)")
+ if(object at family$link=="identity") object at family$linkinv(newx%*%object at parameters$coefficients)
+ else object at family$linkinv(newx%*%object at parameters$coefficients,base=object at family$base)
+ }
+ }
+)
+
+setMethod("fit","transInit",
+ function(object,w,ntimes) {
+ pars <- object at parameters
+ if(missing(w)) w <- NULL
+ pars <- object at parameters
+ base <- object at family$base # delete me
+ #y <- object at y[,-base]
+ y <- object at y
+ x <- object at x
+ if(is.matrix(y)) na <- unlist(apply(y,2,function(x) which(is.na(x)))) else na <- which(is.na(y))
+ if(is.matrix(x)) na <- c(na,unlist(apply(x,2,function(x) which(is.na(x))))) else na <- c(na,which(is.na(x)))
+ if(!is.null(w)) na <- c(na,which(is.na(w)))
+ y <- as.matrix(y)
+ x <- as.matrix(x)
+ na <- unique(na)
+ if(length(na)>0) {
+ x <- x[-na,]
+ y <- y[-na,]
+ #y <- round(y) # delete me
+ if(!is.null(w)) w <- w[-na]
+ }
+ switch(object at family$link,
+ mlogit = {
+ mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
+ mask[,base] <- 0 # fix base category coefficients to 0
+ mask <- rbind(0,mask) # fix "bias" nodes to 0
+ Wts <- mask
+ Wts[-1,] <- pars$coefficients # set starting weights
+ Wts[Wts == Inf] <- .Machine$double.max.exp # Fix this!!!!
+ Wts[Wts == -Inf] <- .Machine$double.min.exp # Fix this!!!!!
+ if(!is.null(w)) {
+ if(NCOL(y) < 3) {
+ fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ } else {
+ fit <- nnet.default(x,y,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ }
+ } else {
+ if(NCOL(y) < 3) {
+ fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ } else {
+ fit <- nnet.default(x,y,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ }
+ }
+ pars$coefficients <- t(matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1)[-1,])
+ object <- setpars(object,unlist(pars))
+ },
+ identity = {
+ # object at y = fbo$xi[,,i]/fbo$gamma[,i]
+ # should become (see em):
+ #for(k in 1:ns) {
+ # trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
+ # }
+ if(!is.null(w)) {
+ sw <- sum(w)
+ pars <- colSums(w*object at y)/sum(w)
+ } else {
+ pars <- colMeans(object at y)
+ }
+ object <- setpars(object,pars)
+ },
+ stop("link function not implemented")
+ )
+ object
+ }
+)
+
+setMethod("simulate",signature(object="transInit"),
+ function(object,nsim=1,seed=NULL,times,is.prior=FALSE,...) {
+ if(!is.null(seed)) set.seed(seed)
+ if(is.prior) {
+ pr <- dens(object)
+ sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nrow(pr)))
+ states <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ return(states)
+ } else {
+ if(missing(times)) {
+ # this is likely to be a stationary model...???
+ pr <- predict(object)
+ } else {
+ pr <- predict(object)[times,]
+ if(length(times)==1) pr <- matrix(pr,ncol=length(pr))
+ }
+ nt <- nrow(pr)
+ sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
+ states <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ # states <- apply(apply(pr,2,rmultinom rmultinom(nt*nsim,size=1,prob=pr),2,function(x) which(x==1))
+ return(states)
+ }
+ }
+)
+
+setMethod("getdf","transInit",
+ function(object) {
+ npar <- sum(!object at fixed)
+ if(object at family$link == "identity") {
+ npar <- npar-1
+ if(npar<0) npar <- 0
+ }
+ npar
+ }
+)
Modified: trunk/depmixNew-test6.R
===================================================================
--- trunk/depmixNew-test6.R 2010-03-08 16:09:22 UTC (rev 393)
+++ trunk/depmixNew-test6.R 2010-03-08 18:47:27 UTC (rev 394)
@@ -56,7 +56,8 @@
allpars <- getpars(fmod1)
-allpars[2]=Inf # this means the process will always start in state 2
+allpars[2]=1 # this means the process will always start in state 2
+allpars[1]=0
allpars[14]=0 # the corr parameters in state 1 are now both 0, corresponding the 0.5 prob
allpars[c(4,8)] <- -4
@@ -69,9 +70,15 @@
conpat[4] <- conpat[8] <- 2
conpat[6] <- conpat[10] <- 3
-fm1sol <- fit(stmod,equal=conpat,method="rsolnp")
+logLik(stmod)
+
+# source("R/depmixfit.R")
+
fm1don <- fit(stmod,equal=conpat,method="donlp")
+
+fm1sol <- fit(stmod,equal=conpat,method="rsolnp")
+
fm1sol
fm1don
Modified: trunk/man/depmix.fit.Rd
===================================================================
--- trunk/man/depmix.fit.Rd 2010-03-08 16:09:22 UTC (rev 393)
+++ trunk/man/depmix.fit.Rd 2010-03-08 18:47:27 UTC (rev 394)
@@ -76,9 +76,10 @@
\details{
Models are fitted by the EM algorithm if there are no constraints on
- the parameters. Otherwise the general optimizers \code{donlp2} (from
- package \code{Rdonlp2}) or \code{solnp} (from package \code{Rsolnp})
- are used which handle general linear (in-)equality constraints.
+ the parameters. Otherwise the general optimizers \code{solnp}, the
+ default (from package \code{Rsolnp}) or \code{donlp2} (from package
+ \code{Rdonlp2}) are used which handle general linear (in-)equality
+ constraints.
Three types of constraints can be specified on the parameters: fixed,
equality, and general linear (in-)equality constraints. Constraint
@@ -91,7 +92,7 @@
estimated to be equal. Any integers can be used in this way except 0
and 1, which indicate fixed and free parameters, respectively.
- Using \code{donlp2} or \code{solnp}, a Newton-Raphson scheme is employed
+ Using \code{solnp} or \code{donlp2} , a Newton-Raphson scheme is employed
to estimate parameters subject to linear constraints by imposing:
bl <= A*x <= bu,
More information about the depmix-commits
mailing list