[Depmix-commits] r391 - in trunk: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 8 16:26:42 CET 2010
Author: ingmarvisser
Date: 2010-03-08 16:26:42 +0100 (Mon, 08 Mar 2010)
New Revision: 391
Modified:
trunk/R/EM.R
trunk/R/allGenerics.R
trunk/R/depmix-class.R
trunk/R/depmixfit.R
trunk/R/freepars.R
trunk/R/response-class.R
trunk/R/responseGLM.R
trunk/R/responseGLMMULTINOM.R
trunk/R/responseMVN.R
trunk/R/transInit.R
trunk/depmixNew-test6.R
trunk/man/GLMresponse.Rd
trunk/man/depmix.fit.Rd
trunk/man/response-class.Rd
trunk/man/transInit.Rd
Log:
Added constraints and parameter bounds to all submodels that have them (sd and multinomial identity models)
Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/EM.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -1,206 +1,202 @@
-#
-# Maarten Speekenbrink 23-3-2008
-#
-
-em <- function(object,maxit=100,tol=1e-8,crit=c(relative,absolute),verbose=FALSE,...) {
- if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
- call <- match.call()
- if(is(object,"depmix")) {
- call[[1]] <- as.name("em.depmix")
- } else {
- call[[1]] <- as.name("em.mix")
- }
- object <- eval(call, parent.frame())
- object
-}
-
-# em for lca and mixture models
-em.mix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),verbose=FALSE,...) {
- if(!is(object,"mix")) stop("object is not of class 'mix'")
+#
+# Maarten Speekenbrink 23-3-2008
+#
+
+em <- function(object,maxit=100,tol=1e-8,crit=c(relative,absolute),verbose=FALSE,...) {
+ if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
+ call <- match.call()
+ if(is(object,"depmix")) {
+ call[[1]] <- as.name("em.depmix")
+ } else {
+ call[[1]] <- as.name("em.mix")
+ }
+ object <- eval(call, parent.frame())
+ object
+}
+
+# em for lca and mixture models
+em.mix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),verbose=FALSE,...) {
+
+ if(!is(object,"mix")) stop("object is not of class 'mix'")
+
crit <- match.arg(crit)
-
- ns <- object at nstates
-
- ntimes <- ntimes(object)
- lt <- length(ntimes)
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
-
- converge <- FALSE
- j <- 0
-
- # compute responsibilities
- B <- apply(object at dens,c(1,3),prod)
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
- # normalize
- gamma <- gamma/rowSums(gamma)
-
- LL.old <- LL + 1
-
- while(j <= maxit & !converge) {
-
- # maximization
-
- # should become object at prior <- fit(object at prior)
- object at prior@y <- gamma[bt,,drop=FALSE]
- object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
- object at init <- dens(object at prior)
-
- for(i in 1:ns) {
- for(k in 1:nresp(object)) {
- object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
- # update dens slot of the model
- object at dens[,k,i] <- dens(object at response[[i]][[k]])
- }
- }
-
- # expectation
- B <- apply(object at dens,c(1,3),prod)
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
- # normalize
- gamma <- gamma/rowSums(gamma)
-
- # print stuff
- if(verbose&((j%%5)==0)) {
- cat("iteration",j,"logLik:",LL,"\n")
- }
-
+
+ ns <- object at nstates
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ converge <- FALSE
+ j <- 0
+
+ # compute responsibilities
+ B <- apply(object at dens,c(1,3),prod)
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+
+ LL.old <- LL + 1
+
+ while(j <= maxit & !converge) {
+
+ # maximization
+
+ # should become object at prior <- fit(object at prior)
+ object at prior@y <- gamma[bt,,drop=FALSE]
+ object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
+ object at init <- dens(object at prior)
+
+ for(i in 1:ns) {
+ for(k in 1:nresp(object)) {
+ object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
+ # update dens slot of the model
+ object at dens[,k,i] <- dens(object at response[[i]][[k]])
+ }
+ }
+
+ # expectation
+ B <- apply(object at dens,c(1,3),prod)
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+
+ # print stuff
+ if(verbose&((j%%5)==0)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ }
+
if(LL >= LL.old) {
- if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
- cat("iteration",j,"logLik:",LL,"\n")
+ if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
converge <- TRUE
- }
+ }
} else {
# this should not really happen...
if(j > 0) warning("likelihood decreased on iteration",j)
- }
-
- LL.old <- LL
- j <- j+1
-
- }
-
- class(object) <- "mix.fitted"
-
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ class(object) <- "mix.fitted"
+
if(converge) {
- object at message <- switch(crit,
- relative = "Log likelihood converged to within tol. (relative change crit.)",
- absolute = "Log likelihood converged to within tol. (absolute change crit.)"
- )
- } else object at message <- "'maxit' iterations reached in EM without convergence."
-
- # no constraints in EM
- object at conMat <- matrix()
- object at lin.lower <- numeric()
- object at lin.upper <- numeric()
-
- object
-
-}
-
-# em for hidden markov models
-em.depmix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),verbose=FALSE,...) {
-
- if(!is(object,"depmix")) stop("object is not of class '(dep)mix'")
+ object at message <- switch(crit,
+ "relative" = "Log likelihood converged to within tol. (relative change crit.)",
+ "absolute" = "Log likelihood converged to within tol. (absolute change crit.)"
+ )
+ } else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM
+ object at conMat <- matrix()
+ object at lin.lower <- numeric()
+ object at lin.upper <- numeric()
+
+ object
+
+}
+
+# em for hidden markov models
+em.depmix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),verbose=FALSE,...) {
+
+ if(!is(object,"depmix")) stop("object is not of class '(dep)mix'")
crit <- match.arg(crit)
-
- ns <- object at nstates
-
- ntimes <- ntimes(object)
- lt <- length(ntimes)
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
-
- converge <- FALSE
- j <- 0
-
- # A <- object at trDens
- # B <- object at dens
- # init <- object at init
-
- # initial expectation
- fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
- LL <- fbo$logLike
- LL.old <- LL + 1
-
- while(j <= maxit & !converge) {
-
- # maximization
-
- # should become object at prior <- fit(object at prior)
- object at prior@y <- fbo$gamma[bt,,drop=FALSE]
- object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
- object at init <- dens(object at prior)
-
- trm <- matrix(0,ns,ns)
- for(i in 1:ns) {
- if(max(ntimes(object)>1)) { # skip transition parameters update in case of latent class model
- if(!object at stationary) {
- object at transition[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
- object at transition[[i]] <- fit(object at transition[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=ntimes(object)) # check this
- } else {
- for(k in 1:ns) {
- trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
- }
- # FIX THIS; it will only work with specific trinModels??
- object at transition[[i]]@parameters$coefficients <- switch(object at transition[[i]]@family$link,
+
+ ns <- object at nstates
+
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ converge <- FALSE
+ j <- 0
+
+ # A <- object at trDens
+ # B <- object at dens
+ # init <- object at init
+
+ # initial expectation
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ LL <- fbo$logLike
+ LL.old <- LL + 1
+
+ while(j <= maxit & !converge) {
+
+ # maximization
+
+ # should become object at prior <- fit(object at prior)
+ object at prior@y <- fbo$gamma[bt,,drop=FALSE]
+ object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
+ object at init <- dens(object at prior)
+
+ trm <- matrix(0,ns,ns)
+ for(i in 1:ns) {
+ if(!object at stationary) {
+ object at transition[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
+ object at transition[[i]] <- fit(object at transition[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=ntimes(object)) # check this
+ } else {
+ for(k in 1:ns) {
+ trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
+ }
+ # FIX THIS; it will only work with specific trinModels??
+ object at transition[[i]]@parameters$coefficients <- switch(object at transition[[i]]@family$link,
identity = object at transition[[i]]@family$linkfun(trm[i,]),
mlogit = object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base),
- object at transition[[i]]@family$linkfun(trm[i,]))
- }
- # update trDens slot of the model
- object at trDens[,,i] <- dens(object at transition[[i]])
- }
- }
-
- for(i in 1:ns) {
-
- for(k in 1:nresp(object)) {
- object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
- # update dens slot of the model
- object at dens[,k,i] <- dens(object at response[[i]][[k]])
- }
- }
-
- # expectation
- fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
- LL <- fbo$logLike
-
+ object at transition[[i]]@family$linkfun(trm[i,])
+ )
+ }
+ # update trDens slot of the model
+ object at trDens[,,i] <- dens(object at transition[[i]])
+ }
+
+ for(i in 1:ns) {
+ for(k in 1:nresp(object)) {
+ object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
+ # update dens slot of the model
+ object at dens[,k,i] <- dens(object at response[[i]][[k]])
+ }
+ }
+
+ # expectation
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ LL <- fbo$logLike
+
if(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
-
+
if( (LL >= LL.old)) {
- if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
- cat("iteration",j,"logLik:",LL,"\n")
+ if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
converge <- TRUE
- }
+ }
} else {
# this should not really happen...
if(j > 0) warning("likelihood decreased on iteration",j)
- }
-
- LL.old <- LL
- j <- j+1
-
- }
-
- #if(class(object)=="depmix") class(object) <- "depmix.fitted"
- #if(class(object)=="mix") class(object) <- "mix.fitted"
-
- class(object) <- "depmix.fitted"
-
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ class(object) <- "depmix.fitted"
+
if(converge) {
- object at message <- switch(crit,
- relative = "Log likelihood converged to within tol. (relative change crit.)",
- absolute = "Log likelihood converged to within tol. (absolute change crit.)"
- )
- } else object at message <- "'maxit' iterations reached in EM without convergence."
-
- # no constraints in EM
- object at conMat <- matrix()
- object at lin.lower <- numeric()
- object at lin.upper <- numeric()
-
- object
-}
+ object at message <- switch(crit,
+ "relative" = "Log likelihood converged to within tol. (relative change crit.)",
+ "absolute" = "Log likelihood converged to within tol. (absolute change crit.)"
+ )
+ } else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM
+ object at conMat <- matrix()
+ object at lin.lower <- numeric()
+ object at lin.upper <- numeric()
+
+ object
+}
Modified: trunk/R/allGenerics.R
===================================================================
--- trunk/R/allGenerics.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/allGenerics.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -8,7 +8,6 @@
require(methods)
require(MASS)
require(nnet)
- require(MCMCpack)
}
.Last.lib <- function(libpath) {}
Modified: trunk/R/depmix-class.R
===================================================================
--- trunk/R/depmix-class.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/depmix-class.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -17,7 +17,7 @@
setClass("mix",
representation(response="list", # response models
- prior="ANY", # the prior model (multinomial logistic)
+ prior="ANY", # the prior model (multinomial)
dens="array", # response densities (B)
init="array", # usually called pi
nstates="numeric",
@@ -65,26 +65,11 @@
# simulate state sequences first, then observations
# random generation is slow when done separately for each t, so first draw
- # variates for all t, and then determine state sequences iteratively
+ # variates for all t, and then determine state sequences iteratively
states <- array(,dim=c(nt,nsim))
states[bt,] <- simulate(object at prior,n=nsim,is.prior=T)
sims <- array(,dim=c(nt,ns,nsim))
-
-# for(i in 1:ns) {
-# if(is.stationary(object)) {
-# # TODO: this is a temporary fix!!!
-# sims[,i,] <- simulate(object at transition[[i]],nsim=nsim,times=rep(1,nt))
-# } else {
-# sims[,i,] <- simulate(object at transition[[i]],nsim=nsim)
-# }
-# }
-# # track states
-# for(case in 1:lt) {
-# for(i in (bt[case]+1):et[case]) {
-# states[i,] <- sims[cbind(i,states[i-1,],1:nsim)]
-# }
-# }
-
+
states <- as.vector(states)
responses <- list(length=nr)
#responses <- array(,dim=c(nt,nr,nsim))
@@ -102,7 +87,6 @@
object at prior@x <- as.matrix(apply(object at prior@x,2,rep,nsim))
for(j in 1:ns) {
-# if(!is.stationary(object)) object at transition[[j]]@x <- as.matrix(apply(object at transition[[j]]@x,2,rep,nsim))
for(i in 1:nr) {
object at response[[j]][[i]]@y <- as.matrix(responses[[i]])
object at response[[j]][[i]]@x <- as.matrix(apply(object at response[[j]][[i]]@x,2,rep,nsim))
Modified: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/depmixfit.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -13,7 +13,7 @@
# otherwise EM is good
if(is.null(method)) {
if(constr) {
- method="donlp"
+ method="rsolnp"
} else {
method="EM"
}
@@ -57,28 +57,77 @@
# get the full set of parameters
allpars <- getpars(object)
- # get the reduced set of parameters, ie the ones that will be optimized
+
+ # get the reduced set of parameters, ie the ones that will be optimized
pars <- allpars[!fixed]
# set bounds, if any (should add bounds for eg sd parameters at some point ...)
- par.u <- rep(+Inf, length(pars))
- par.l <- rep(-Inf, length(pars))
+ par.u <- rep(+Inf, npar(object))
+ par.l <- rep(-Inf, npar(object))
- # make loglike function that only depends on pars
- logl <- function(pars) {
- allpars[!fixed] <- pars
- object <- setpars(object,allpars)
- ans = -as.numeric(logLik(object))
- if(is.na(ans)) ans = 100000
- ans
- }
-
# make constraint matrix and its upper and lower bounds
lincon <- matrix(0,nr=0,nc=npar(object))
lin.u <- numeric(0)
lin.l <- numeric(0)
- # incorporate equality constraints, if any
+ ns <- nstates(object)
+ nrsp <- nresp(object)
+
+ # get bounds from submodels
+ # get internal linear constraints from submodels
+
+ # first for the prior model
+ bp <- 1
+ ep <- npar(object at prior)
+ if(!is.null(object at prior@constr)) {
+ par.u[bp:ep] <- object at prior@constr$parup
+ par.l[bp:ep] <- object at prior@constr$parlow
+ # add linear constraints, if any
+ if(!is.null(object at prior@constr$lin)) {
+ lincon <- rbind(lincon,0)
+ lincon[nrow(lincon),bp:ep] <- object at prior@constr$lin
+ lin.u[nrow(lincon)] <- object at prior@constr$linup
+ lin.l[nrow(lincon)] <- object at prior@constr$linlow
+ }
+ }
+
+ # ... for the transition models
+ if(is(object,"depmix")) {
+ for(i in 1:ns) {
+ bp <- ep + 1
+ ep <- ep+npar(object at transition[[i]])
+ if(!is.null(object at transition[[i]]@constr)) {
+ par.u[bp:ep] <- object at transition[[i]]@constr$parup
+ par.l[bp:ep] <- object at transition[[i]]@constr$parlow
+ }
+ if(!is.null(object at transition[[i]]@constr$lin)) {
+ lincon <- rbind(lincon,0)
+ lincon[nrow(lincon),bp:ep] <- object at transition[[i]]@constr$lin
+ lin.u[nrow(lincon)] <- object at transition[[i]]@constr$linup
+ lin.l[nrow(lincon)] <- object at transition[[i]]@constr$linlow
+ }
+ }
+ }
+
+ # ... for the response models
+ for(i in 1:ns) {
+ for(j in 1:nrsp) {
+ bp <- ep + 1
+ ep <- ep + npar(object at response[[i]][[j]])
+ if(!is.null(object at response[[i]][[j]]@constr)) {
+ par.u[bp:ep] <- object at response[[i]][[j]]@constr$parup
+ par.l[bp:ep] <- object at response[[i]][[j]]@constr$parlow
+ }
+ if(!is.null(object at response[[i]][[j]]@constr$lin)) {
+ lincon <- rbind(lincon,0)
+ lincon[nrow(lincon),bp:ep] <- object at response[[i]][[j]]@constr$lin
+ lin.u[nrow(lincon)] <- object at response[[i]][[j]]@constr$linup
+ lin.l[nrow(lincon)] <- object at response[[i]][[j]]@constr$linlow
+ }
+ }
+ }
+
+ # incorporate equality constraints provided with the fit function, if any
if(eq) {
if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
equal <- pa2conr(equal)$conr
@@ -105,10 +154,43 @@
}
}
+ print(round(allpars,2))
+ print(par.u)
+ print(par.l)
+
+ print(lincon)
+ print(lin.u)
+ print(lin.l)
+
# select only those columns of the constraint matrix that correspond to non-fixed parameters
linconFull <- lincon
- lincon <- lincon[,!fixed,drop=FALSE]
+ lincon <- lincon[,!fixed,drop=FALSE]
+
+ # remove redundant rows in lincon (all zeroes)
+ allzero <- which(apply(lincon,1,function(y) all(y==0)))
+ if(length(allzero)>0) {
+ lincon <- lincon[-allzero,,drop=FALSE]
+ lin.u <- lin.u[-allzero]
+ lin.l <- lin.l[-allzero]
+ }
+ print(round(pars,2))
+ print(par.u[!fixed])
+ print(par.l[!fixed])
+
+ print(lincon)
+ print(lin.u)
+ print(lin.l)
+
+ # make loglike function that only depends on pars
+ logl <- function(pars) {
+ allpars[!fixed] <- pars
+ object <- setpars(object,allpars)
+ ans = -as.numeric(logLik(object))
+ if(is.na(ans)) ans = 100000
+ ans
+ }
+
if(method=="donlp") {
# set donlp2 control parameters
cntrl <- donlp2.control(hessian=FALSE,difftype=2,report=TRUE)
@@ -119,8 +201,8 @@
# optimize the parameters
result <- donlp2(pars,logl,
- par.upper=par.u,
- par.lower=par.l,
+ par.upper=par.u[!fixed],
+ par.lower=par.l[!fixed],
A=lincon,
lin.upper=lin.u,
lin.lower=lin.l,
@@ -146,7 +228,9 @@
ineq <- which(lin.u!=lin.l)
if(length(ineq)>0) lineq <- lincon[-ineq, ,drop=FALSE]
else lineq <- lincon
-
+
+ print(lincon)
+
# returns the evaluated equality constraints
if(nrow(lineq)>0) {
eqfun <- function(pp) {
@@ -161,6 +245,8 @@
lineq.bound=NULL
}
+ print(lineq.bound)
+
# select the inequality constraints
if(length(ineq)>0) {
linineq <- lincon[ineq, ,drop=FALSE]
@@ -182,8 +268,8 @@
ineqLB = ineqLB,
ineqUB = ineqUB,
ineqgrad = NULL,
- LB = NULL,
- UB = NULL,
+ LB = par.l[!fixed],
+ UB = par.u[!fixed],
control = list(trace = 1)
)
Modified: trunk/R/freepars.R
===================================================================
--- trunk/R/freepars.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/freepars.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -1,11 +1,34 @@
-# depends on getpars(object)
+# depends on getpars(object) & getdf
setMethod("freepars","mix",
function(object) {
- free <- sum(!getpars(object,which="fixed"))
+ free <- getdf(object at prior)
+ ns <- nstates(object)
+ nresp <- nresp(object)
+ for(i in 1:ns) {
+ for(j in 1: nresp) {
+ free <- free + getdf(object at response[[i]][[j]])
+ }
+ }
free
}
)
+# depends on getpars(object) & getdf
+setMethod("freepars","depmix",
+ function(object) {
+ free <- getdf(object at prior)
+ ns <- nstates(object)
+ nresp <- nresp(object)
+ for(i in 1:ns) {
+ free <- free + getdf(object at transition[[i]])
+ for(j in 1: nresp) {
+ free <- free + getdf(object at response[[i]][[j]])
+ }
+ }
+ free
+ }
+)
+
# depends on nlin(object) and getpars(object)
setMethod("freepars","depmix.fitted",
function(object) {
Modified: trunk/R/response-class.R
===================================================================
--- trunk/R/response-class.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/response-class.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -35,5 +35,3 @@
}
)
-
-
Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/responseGLM.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -32,6 +32,10 @@
parameters$coefficients <- vector("numeric",length=ncol(x))
if(family$family=="gaussian") {
parameters$sd <- 1
+ constr <- list(
+ parup = rep(Inf,ncol(y)+1),
+ parlow = c(rep(-Inf,ncol(y)),0)
+ )
}
if(family$family=="binomial") {
# FIX ME
Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/responseGLMMULTINOM.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -85,3 +85,14 @@
return(sims)
}
)
+
+setMethod("getdf","MULTINOMresponse",
+ function(object) {
+ npar <- sum(!object at fixed)
+ if(object at family$link == "identity") {
+ npar <- npar-1
+ if(npar<0) npar <- 0
+ }
+ npar
+ }
+)
Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/responseMVN.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -16,187 +16,188 @@
}
-
-setClass("MVNresponse",
- representation(formula="formula"),
- contains="response"
-)
-
-setMethod("fit","MVNresponse",
- function(object,w) {
- if(missing(w)) w <- NULL
- pars <- object at parameters
- if(!is.null(w)) fit <- lm.wfit(x=object at x,y=object at y,w=w) else fit <- lm.fit(x=object at x,y=object at y)
- object at parameters$coefficients <- fit$coefficients
- if(!is.null(w)) object at parameters$Sigma <- cov2par(cov.wt(x=fit$residuals,wt=w)$cov) else object at parameters$Sigma <- cov2par(cov(fit$residuals))
- object
- }
-)
-
-dm_dmvnorm <- function(x,mean,sigma,log=FALSE,logdet,invSigma) {
- # taken from mvtnorm package
- # allows passing of logdet (sigma) and invsigma to save
- # computation when called repeated times with same sigma
- if (is.vector(x)) {
- x <- matrix(x, ncol = length(x))
- }
- if (missing(mean)) {
- mean <- rep(0, length = ncol(x))
- }
- if(missing(invSigma)) {
- if (missing(sigma)) {
- sigma <- diag(ncol(x))
- }
- invSigma <- solve(sigma)
- }
- # check consistency
-
- if (NCOL(x) != NCOL(invSigma)) {
- stop("x and sigma have non-conforming size")
- }
- if (NROW(invSigma) != NCOL(invSigma)) {
- stop("sigma must be a square matrix")
- }
- if (NCOL(invSigma) != NCOL(mean)) {
- stop("mean and sigma have non-conforming size")
- }
- if(NROW(mean) == NROW(x)) {
- # varying means
-
- # from "mahalanobis":
- x <- as.matrix(x) - mean
- distval <- rowSums((x %*% invSigma) * x)
- #names(retval) <- rownames(x)
- #retval
- } else {
- # constant mean
- if (length(mean) != NROW(invSigma)) {
- stop("mean and sigma have non-conforming size")
- }
- distval <- mahalanobis(x, center = mean, cov = invSigma, inverted=TRUE)
- }
- if(missing(logdet)) logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
- logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
- if (log) {
- return(logretval)
- } else {
- return(exp(logretval))
- }
-}
-
-
-setMethod("logDens","MVNresponse",
- function(object,...) {
- dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=TRUE,...)
- }
-)
-
-setMethod("dens","MVNresponse",
- function(object,log=FALSE,...) {
- dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=log,...)
- }
-)
-
-
-setMethod("predict","MVNresponse",
- function(object) {
- object at x%*%object at parameters$coefficients
- }
-)
-
-setMethod("simulate",signature(object="MVNresponse"),
- function(object,nsim=1,seed=NULL,times) {
- if(!is.null(seed)) set.seed(seed)
- if(missing(times)) {
- # draw in one go
- mu <- predict(object)
- } else {
- mu <- predict(object)[times,]
- }
- nt <- nrow(mu)
- if(nrow(object at parameters$coefficients==1)) response <- mvrnorm(nt*nsim,mu=mu[1,],Sigma=par2cov(object at parameters$Sigma))
- else {
- response <- matrix(0,nrow(mu),ncol(mu))
- for(i in 1:nrow(mu)) {
- response[i,] <- response <- mvrnorm(1,mu=mu[i,],Sigma=par2cov(object at parameters$Sigma))
- }
- }
- return(response)
- }
-)
-
-setMethod("MVNresponse",
- signature(formula="formula"),
- function(formula,data,pstart=NULL,fixed=NULL,...) {
- call <- match.call()
- 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 <- model.response(mf)
- if(!is.matrix(y)) y <- matrix(y,ncol=1)
- parameters <- list()
- parameters$coefficients <- matrix(0.0,ncol=ncol(y),nrow=ncol(x))
- parameters$Sigma <- cov2par(diag(ncol(y)))
- npar <- length(unlist(parameters))
- if(is.null(fixed)) fixed <- as.logical(rep(0,npar))
- if(!is.null(pstart)) {
- if(length(pstart)!=npar) stop("length of 'pstart' must be",npar)
- parameters$coefficients <- matrix(pstart[1:(ncol(x)*ncol(y))],ncol(x),byrow=T)
- parameters$Sigma <- as.numeric(pstart[(length(parameters$coefficients)+1):length(pstart)])
- }
- mod <- new("MVNresponse",formula=formula,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
- mod
- }
-)
-
-setMethod("show","MVNresponse",
- function(object) {
- cat("Multivariate Normal Model, formula: ",sep="")
- print(object at formula)
- cat("Coefficients: \n")
- print(object at parameters$coefficients)
- cat("Sigma: \n")
- print(par2cov(object at parameters$Sigma))
- }
-)
-
-setMethod("setpars","MVNresponse",
- function(object, values, which="pars", prob=FALSE, ...) {
- npar <- npar(object)
- if(length(values)!=npar) stop("length of 'values' must be",npar)
- # determine whether parameters or fixed constraints are being set
- nms <- names(object at parameters$coefficients)
- switch(which,
- "pars" = {
- object at parameters$coefficients <- matrix(values[1:length(object at parameters$coefficients)],ncol(object at x))
- st <- length(object at parameters$coefficients)+1
- object at parameters$Sigma <- as.numeric(values[st:(st+length(object at parameters$Sigma)-1)])
- },
- "fixed" = {
- object at fixed <- as.logical(values)
- }
- )
- names(object at parameters$coefficients) <- nms
- return(object)
- }
-)
-
-setMethod("getpars","MVNresponse",
- function(object,which="pars",...) {
- switch(which,
- "pars" = {
- parameters <- numeric()
- parameters <- unlist(object at parameters)
- pars <- parameters
- },
- "fixed" = {
- pars <- object at fixed
- }
- )
- return(pars)
- }
-)
+
+setClass("MVNresponse",
+ representation(formula="formula"),
+ contains="response"
+)
+
+setMethod("fit","MVNresponse",
+ function(object,w) {
+ if(missing(w)) w <- NULL
+ pars <- object at parameters
+ if(!is.null(w)) fit <- lm.wfit(x=object at x,y=object at y,w=w) else fit <- lm.fit(x=object at x,y=object at y)
+ object at parameters$coefficients <- fit$coefficients
+ if(!is.null(w)) object at parameters$Sigma <- cov2par(cov.wt(x=fit$residuals,wt=w)$cov) else object at parameters$Sigma <- cov2par(cov(fit$residuals))
+ object
+ }
+)
+
+dm_dmvnorm <- function(x,mean,sigma,log=FALSE,logdet,invSigma) {
+ # taken from mvtnorm package
+ # allows passing of logdet (sigma) and invsigma to save
+ # computation when called repeated times with same sigma
+ if (is.vector(x)) {
+ x <- matrix(x, ncol = length(x))
+ }
+ if (missing(mean)) {
+ mean <- rep(0, length = ncol(x))
+ }
+ if(missing(invSigma)) {
+ if (missing(sigma)) {
+ sigma <- diag(ncol(x))
+ }
+ invSigma <- solve(sigma)
+ }
+ # check consistency
+
+ if (NCOL(x) != NCOL(invSigma)) {
+ stop("x and sigma have non-conforming size")
+ }
+ if (NROW(invSigma) != NCOL(invSigma)) {
+ stop("sigma must be a square matrix")
+ }
+ if (NCOL(invSigma) != NCOL(mean)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ if(NROW(mean) == NROW(x)) {
+ # varying means
+
+ # from "mahalanobis":
+ x <- as.matrix(x) - mean
+ distval <- rowSums((x %*% invSigma) * x)
+ #names(retval) <- rownames(x)
+ #retval
+ } else {
+ # constant mean
+ if (length(mean) != NROW(invSigma)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ distval <- mahalanobis(x, center = mean, cov = invSigma, inverted=TRUE)
+ }
+ if(missing(logdet)) logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
+ logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
+ if (log) {
+ return(logretval)
+ } else {
+ return(exp(logretval))
+ }
+}
+
+
+setMethod("logDens","MVNresponse",
+ function(object,...) {
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=TRUE,...)
+ }
+)
+
+setMethod("dens","MVNresponse",
+ function(object,log=FALSE,...) {
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=log,...)
+ }
+)
+
+
+setMethod("predict","MVNresponse",
+ function(object) {
+ object at x%*%object at parameters$coefficients
+ }
+)
+
+setMethod("simulate",signature(object="MVNresponse"),
+ function(object,nsim=1,seed=NULL,times) {
+ if(!is.null(seed)) set.seed(seed)
+ if(missing(times)) {
+ # draw in one go
+ mu <- predict(object)
+ } else {
+ mu <- predict(object)[times,]
+ }
+ nt <- nrow(mu)
+ if(nrow(object at parameters$coefficients==1)) response <- mvrnorm(nt*nsim,mu=mu[1,],Sigma=par2cov(object at parameters$Sigma))
+ else {
+ response <- matrix(0,nrow(mu),ncol(mu))
+ for(i in 1:nrow(mu)) {
+ response[i,] <- response <- mvrnorm(1,mu=mu[i,],Sigma=par2cov(object at parameters$Sigma))
+ }
+ }
+ return(response)
+ }
+)
+
+setMethod("MVNresponse",
+ signature(formula="formula"),
+ function(formula,data,pstart=NULL,fixed=NULL,...) {
+ call <- match.call()
+ 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 <- model.response(mf)
+ if(!is.matrix(y)) y <- matrix(y,ncol=1)
+ parameters <- list()
+ constr <- NULL
+ parameters$coefficients <- matrix(0.0,ncol=ncol(y),nrow=ncol(x))
+ parameters$Sigma <- cov2par(diag(ncol(y)))
+ npar <- length(unlist(parameters))
+ if(is.null(fixed)) fixed <- as.logical(rep(0,npar))
+ if(!is.null(pstart)) {
+ if(length(pstart)!=npar) stop("length of 'pstart' must be",npar)
+ parameters$coefficients <- matrix(pstart[1:(ncol(x)*ncol(y))],ncol(x),byrow=T)
+ parameters$Sigma <- as.numeric(pstart[(length(parameters$coefficients)+1):length(pstart)])
+ }
+ mod <- new("MVNresponse",formula=formula,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr)
+ mod
+ }
+)
+
+setMethod("show","MVNresponse",
+ function(object) {
+ cat("Multivariate Normal Model, formula: ",sep="")
+ print(object at formula)
+ cat("Coefficients: \n")
+ print(object at parameters$coefficients)
+ cat("Sigma: \n")
+ print(par2cov(object at parameters$Sigma))
+ }
+)
+
+setMethod("setpars","MVNresponse",
+ function(object, values, which="pars", prob=FALSE, ...) {
+ npar <- npar(object)
+ if(length(values)!=npar) stop("length of 'values' must be",npar)
+ # determine whether parameters or fixed constraints are being set
+ nms <- names(object at parameters$coefficients)
+ switch(which,
+ "pars" = {
+ object at parameters$coefficients <- matrix(values[1:length(object at parameters$coefficients)],ncol(object at x))
+ st <- length(object at parameters$coefficients)+1
+ object at parameters$Sigma <- as.numeric(values[st:(st+length(object at parameters$Sigma)-1)])
+ },
+ "fixed" = {
+ object at fixed <- as.logical(values)
+ }
+ )
+ names(object at parameters$coefficients) <- nms
+ return(object)
+ }
+)
+
+setMethod("getpars","MVNresponse",
+ function(object,which="pars",...) {
+ switch(which,
+ "pars" = {
+ parameters <- numeric()
+ parameters <- unlist(object at parameters)
+ pars <- parameters
+ },
+ "fixed" = {
+ pars <- object at fixed
+ }
+ )
+ return(pars)
+ }
+)
Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R 2010-03-05 15:56:16 UTC (rev 390)
+++ trunk/R/transInit.R 2010-03-08 15:26:42 UTC (rev 391)
@@ -23,15 +23,23 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 391
More information about the depmix-commits
mailing list