[Depmix-commits] r674 - in pkg/depmixS4: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 19 17:25:04 CEST 2018
Author: ingmarvisser
Date: 2018-10-19 17:25:03 +0200 (Fri, 19 Oct 2018)
New Revision: 674
Added:
pkg/depmixS4/R/vcov.R
Modified:
pkg/depmixS4/NAMESPACE
pkg/depmixS4/R/hessian.R
Log:
=added vcov function, moved lincon code from hessian to vcov
Modified: pkg/depmixS4/NAMESPACE
===================================================================
--- pkg/depmixS4/NAMESPACE 2018-10-17 15:50:42 UTC (rev 673)
+++ pkg/depmixS4/NAMESPACE 2018-10-19 15:25:03 UTC (rev 674)
@@ -67,7 +67,8 @@
summary,
logLik,
getmodel,
- hessian
+ hessian,
+ vcov
)
useDynLib(depmixS4, .registration = TRUE)
\ No newline at end of file
Modified: pkg/depmixS4/R/hessian.R
===================================================================
--- pkg/depmixS4/R/hessian.R 2018-10-17 15:50:42 UTC (rev 673)
+++ pkg/depmixS4/R/hessian.R 2018-10-19 15:25:03 UTC (rev 674)
@@ -25,81 +25,27 @@
# elements: vector of length npar(object) indicating for each parameter whether it
# is 'inc'luded, 'fix'ed, or estimated on the boundary, 'bnd'; the dimension of the hessian
# is thus the number of non-fixed parameters minus the number of boundary parameters.
-#
-# lincon: the linear constraint matrix needed to compute the variance-covariance
-# matrix; it only contains the parts of the linear constraint matrix that relate to equality
-# constraints; moreover, the columns related to 'fix'ed and boundary ('bnd') parameters
-# are left out.
#
setMethod("hessian", "mix",
- function(object, fixed=NULL, equal=NULL,
- conrows=NULL, conrows.upper=NULL, conrows.lower=NULL,
+ function(object,
tolerance=1e-9,
method="finiteDifferences", ...) {
if(is.nan(logLik(object))) stop("Log likelihood is 'NaN'; cannot compute hessian. ")
+
+ # get the full set of parameters
+ allpars <- getpars(object)
- # check for presence of constraints
- fi <- !is.null(fixed)
- cr <- !is.null(conrows)
- eq <- !is.null(equal)
-
- constr <- any(c(fi,cr,eq))
-
- # determine which parameters are fixed
- if(fi) {
- if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
- } else {
- if(eq) {
- if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
- fixed <- !pa2conr(equal)$free
- } else {
- fixed <- getpars(object,"fixed")
- }
- }
-
- # set those fixed parameters in the appropriate submodels
- object <- setpars(object,fixed,which="fixed")
-
+ # get parameter boundaries from the model object
constraints <- getConstraints(object)
- 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 hessian function, if any
- if(eq) {
- if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
- equal <- pa2conr(equal)$conr
- lincon <- rbind(lincon,equal)
- lin.u <- c(lin.u,rep(0,nrow(equal)))
- lin.l <- c(lin.l,rep(0,nrow(equal)))
- }
+ # get fixed parameters
+ fixed <- getpars(object,"fixed")
- # incorporate general linear constraints, if any, in lincon matrix
- if(cr) {
- if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
- lincon <- rbind(lincon,conrows)
- if(any(conrows.upper==0)) {
- lin.u <- c(lin.u,rep(0,nrow(conrows)))
- } else {
- if(length(conrows.upper)!=nrow(conrows)) stop("'conrows.upper does not have correct length")
- lin.u <- c(lin.u,conrows.upper)
- }
- if(any(conrows.lower==0)) {
- lin.l <- c(lin.l,rep(0,nrow(conrows)))
- } else {
- if(length(conrows.lower)!=nrow(conrows)) stop("'conrows.lower does not have correct length")
- lin.l <- c(lin.l,conrows.lower)
- }
- }
-
- # get the full set of parameters
- allpars <- getpars(object)
-
# return vector with specification of 'inc'luded, 'fix'ed and 'bnd'ary parameters
elements <- rep("inc",npar(object))
@@ -108,34 +54,14 @@
up <- which(sapply(as.numeric(allpars-par.u),all.equal,tolerance=tolerance,0)==TRUE)
bnd <- union(low, up)
-
+ # identify parameters that are fixed or on the boundary
if(length(which(fixed)>0)) elements[which(fixed)] <- "fix"
if(length(bnd)>0) elements[bnd] <- "bnd"
# get the reduced set of parameters, ie the ones that the hessian will be computed for
# only non-fixed parameters
- pars <- allpars[which(elements=="inc")]
-
- # select only those columns of the constraint matrix that correspond to non-fixed parameters
- lincon <- lincon[,which(elements=="inc"),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]
- }
-
- # remove rows of lincon with inequality constraints
- dflu <- lin.u-lin.l
- ineq <- which(dflu!=0)
- if(length(ineq)>0) {
- lincon <- lincon[-ineq,,drop=FALSE]
- lin.u <- lin.u[-ineq]
- lin.l <- lin.l[-ineq]
- }
-
+ pars <- allpars[which(elements=="inc")]
+
# make loglike function that only depends on pars
logl <- function(pars) {
allpars[which(elements=="inc")] <- pars
@@ -146,11 +72,7 @@
}
fdh <- fdHess(pars,logl)
-
- # also return list of length npar that specifies for which parameters
- # the hessian has been computed and for which this has been skipped due
- # to being 1) fixed or 2) on the boundary
- return(list(hessian=fdh$Hessian,elements=elements,lincon=lincon))
+ return(list(hessian=fdh$Hessian,elements=elements))
}
)
\ No newline at end of file
Added: pkg/depmixS4/R/vcov.R
===================================================================
--- pkg/depmixS4/R/vcov.R (rev 0)
+++ pkg/depmixS4/R/vcov.R 2018-10-19 15:25:03 UTC (rev 674)
@@ -0,0 +1,144 @@
+#
+# Ingmar Visser, 19-10-2018
+#
+# Description
+#
+# vcov computes the variance-covariance matrix of depmix model parameters
+#
+# Details
+#
+# The variance-covariance matrix of the parameters is computed using the hessian
+# and the linear constraints in the model. First, the hessian is computed through
+# the 'hessian' method for (dep)mix model. Second, the hessian and linear constraint
+# matrix are passed to the 'hessian2vcov' function.
+#
+# Value
+#
+# The function returns a named list with the following elements:
+#
+# vcov: a square matrix of dimension the 'inc'luded parameters.
+#
+# elements: a vector of length npar(object) indicating
+# for each parameter whether it is 'inc'luded, 'fix'ed, or estimated on the
+# boundary, 'bnd'; the dimension of the vcov matrix is thus the number of
+# non-fixed parameters minus the number of boundary parameters.
+#
+# lincon: the linear constraint matrix used to compute the variance-covariance
+# matrix from the hessian; it only contains the parts of the linear constraint
+# matrix that relate to equality constraints; moreover, the columns related
+# to 'fix'ed and boundary ('bnd') parameters are left out.
+#
+
+setMethod("vcov", "mix",
+ function(object, fixed=NULL, equal=NULL,
+ conrows=NULL, conrows.upper=NULL, conrows.lower=NULL,
+ tolerance=1e-9,
+ method="finiteDifferences", ...) {
+
+ if(is.nan(logLik(object))) stop("Log likelihood is 'NaN'; cannot compute variance-covariance. ")
+
+ # check for presence of constraints
+ fi <- !is.null(fixed)
+ cr <- !is.null(conrows)
+ eq <- !is.null(equal)
+
+ constr <- any(c(fi,cr,eq))
+
+ # determine which parameters are fixed
+ if(fi) {
+ if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
+ } else {
+ if(eq) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ fixed <- !pa2conr(equal)$free
+ } else {
+ fixed <- getpars(object,"fixed")
+ }
+ }
+
+ # set those fixed parameters in the appropriate submodels
+ object <- setpars(object,fixed,which="fixed")
+
+ # get standard constraints from (sub)models
+ constraints <- getConstraints(object)
+
+ 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 hessian function, if any
+ if(eq) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ equal <- pa2conr(equal)$conr
+ lincon <- rbind(lincon,equal)
+ lin.u <- c(lin.u,rep(0,nrow(equal)))
+ lin.l <- c(lin.l,rep(0,nrow(equal)))
+ }
+
+ # incorporate general linear constraints, if any, via argument conrows
+ if(cr) {
+ if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
+ lincon <- rbind(lincon,conrows)
+ if(any(conrows.upper==0)) {
+ lin.u <- c(lin.u,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.upper)!=nrow(conrows)) stop("'conrows.upper does not have correct length")
+ lin.u <- c(lin.u,conrows.upper)
+ }
+ if(any(conrows.lower==0)) {
+ lin.l <- c(lin.l,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.lower)!=nrow(conrows)) stop("'conrows.lower does not have correct length")
+ lin.l <- c(lin.l,conrows.lower)
+ }
+ }
+
+ # get the full set of parameters
+ allpars <- getpars(object)
+
+ # return vector with specification of 'inc'luded, 'fix'ed and 'bnd'ary parameters
+ elements <- rep("inc",npar(object))
+
+ # identify parameters that are on their boundary
+ low <- which(sapply(as.numeric(allpars-par.l),all.equal,tolerance=tolerance,0)==TRUE)
+ up <- which(sapply(as.numeric(allpars-par.u),all.equal,tolerance=tolerance,0)==TRUE)
+ bnd <- union(low, up)
+
+ # identify parameters that are fixed or on the boundary
+ if(length(which(fixed)>0)) elements[which(fixed)] <- "fix"
+ if(length(bnd)>0) elements[bnd] <- "bnd"
+
+ # get the reduced set of parameters, ie the ones that the hessian will be computed for
+ # only non-fixed parameters
+ pars <- allpars[which(elements=="inc")]
+
+ # select only those columns of the constraint matrix that correspond to non-fixed parameters
+ lincon <- lincon[,which(elements=="inc"),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]
+ }
+
+ # remove rows of lincon with inequality constraints
+ dflu <- lin.u-lin.l
+ ineq <- which(dflu!=0)
+ if(length(ineq)>0) {
+ lincon <- lincon[-ineq,,drop=FALSE]
+ lin.u <- lin.u[-ineq]
+ lin.l <- lin.l[-ineq]
+ }
+
+ hs <- hessian(object)
+
+ vc <- hessian2vcov(hs$hessian,lincon)
+
+ return(list(vcov=vc,elements=elements,lincon=lincon))
+
+}
+)
\ No newline at end of file
More information about the depmix-commits
mailing list