[Depmix-commits] r393 - in trunk: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 8 17:09:22 CET 2010
Author: ingmarvisser
Date: 2010-03-08 17:09:22 +0100 (Mon, 08 Mar 2010)
New Revision: 393
Modified:
trunk/NAMESPACE
trunk/R/EM.R
trunk/R/allGenerics.R
trunk/R/depmixfit.R
Log:
More updates for getting getdf right, not quite there yet ...
Modified: trunk/NAMESPACE
===================================================================
--- trunk/NAMESPACE 2010-03-08 15:56:02 UTC (rev 392)
+++ trunk/NAMESPACE 2010-03-08 16:09:22 UTC (rev 393)
@@ -32,7 +32,8 @@
exportMethods(
AIC,
BIC,
- fit,
+ fit,
+ getConstraints,
npar,
freepars,
nlin,
Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R 2010-03-08 15:56:02 UTC (rev 392)
+++ trunk/R/EM.R 2010-03-08 16:09:22 UTC (rev 393)
@@ -92,10 +92,12 @@
)
} 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()
+ # no constraints in EM, except for the standard constraints ...
+ # which are produced by the following (only necessary for getting df right in logLik and such)
+ constraints <- getConstraints(object)
+ object at conMat <- constraints$lincon
+ object at lin.lower <- constraints$lin.l
+ object at lin.upper <- constraints$lin.u
object
@@ -193,10 +195,12 @@
)
} 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()
+ # no constraints in EM, except for the standard constraints ...
+ # which are produced by the following (only necessary for getting df right in logLik and such)
+ constraints <- getConstraints(object)
+ object at conMat <- constraints$lincon
+ object at lin.lower <- constraints$lin.l
+ object at lin.upper <- constraints$lin.u
object
}
Modified: trunk/R/allGenerics.R
===================================================================
--- trunk/R/allGenerics.R 2010-03-08 15:56:02 UTC (rev 392)
+++ trunk/R/allGenerics.R 2010-03-08 16:09:22 UTC (rev 393)
@@ -45,6 +45,8 @@
setGeneric("fit", function(object, ...) standardGeneric("fit"))
+setGeneric("getConstraints", function(object, ...) standardGeneric("getConstraints"))
+
setGeneric("posterior", function(object, ...) standardGeneric("posterior"))
setGeneric("forwardbackward", function(object, ...) standardGeneric("forwardbackward"))
Modified: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2010-03-08 15:56:02 UTC (rev 392)
+++ trunk/R/depmixfit.R 2010-03-08 16:09:22 UTC (rev 393)
@@ -61,71 +61,13 @@
# 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, npar(object))
- par.l <- rep(-Inf, npar(object))
-
- # 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)
+ constraints <- getConstraints(object)
- 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
- }
- }
- }
+ lincon=constraints$lincon
+ lin.u=constraints$lin.u
+ lin.l=constraints$lin.l
+ par.u=constraints$par.u
+ par.l=constraints$par.l
# incorporate equality constraints provided with the fit function, if any
if(eq) {
@@ -154,14 +96,6 @@
}
}
- 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]
@@ -174,14 +108,6 @@
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
@@ -263,11 +189,9 @@
logl,
eqfun = eqfun,
eqB = lineq.bound,
- eqgrad =NULL,
ineqfun = ineqfun,
ineqLB = ineqLB,
ineqUB = ineqUB,
- ineqgrad = NULL,
LB = par.l[!fixed],
UB = par.u[!fixed],
control = list(trace = 1)
@@ -295,3 +219,80 @@
return(object)
}
)
+
+
+setMethod("getConstraints",
+ signature(object="mix"),
+ function(object) {
+
+ # set bounds, if any (should add bounds for eg sd parameters at some point ...)
+ par.u <- rep(+Inf, npar(object))
+ par.l <- rep(-Inf, npar(object))
+
+ # 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)
+
+ 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
+ }
+ }
+ }
+
+ res <- list(lincon=lincon,lin.u=lin.u,lin.l=lin.l,par.u=par.u,par.l=par.l)
+ res
+
+ }
+)
More information about the depmix-commits
mailing list