[Depmix-commits] r670 - in pkg/depmixS4: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 17 11:53:34 CEST 2018
Author: ingmarvisser
Date: 2018-10-17 11:53:34 +0200 (Wed, 17 Oct 2018)
New Revision: 670
Added:
pkg/depmixS4/R/hessian2vcov.R
Removed:
pkg/depmixS4/R/hess2cov.R
Modified:
pkg/depmixS4/DESCRIPTION
pkg/depmixS4/NAMESPACE
pkg/depmixS4/R/allGenerics.R
pkg/depmixS4/R/fdHessHmm.R
Log:
=changed filename from hess2cov.R to hessian2vcov.R
Modified: pkg/depmixS4/DESCRIPTION
===================================================================
--- pkg/depmixS4/DESCRIPTION 2018-09-20 12:59:46 UTC (rev 669)
+++ pkg/depmixS4/DESCRIPTION 2018-10-17 09:53:34 UTC (rev 670)
@@ -4,7 +4,7 @@
Title: Dependent Mixture Models - Hidden Markov Models of GLMs and Other Distributions in S4
Author: Ingmar Visser <i.visser at uva.nl>, Maarten Speekenbrink <m.speekenbrink at ucl.ac.uk>
Maintainer: Ingmar Visser <i.visser at uva.nl>
-Depends: R (>= 3.5.0), nnet, MASS, Rsolnp
+Depends: R (>= 3.5.0), nnet, MASS, Rsolnp, nlme
Imports: stats, stats4, methods
Suggests: gamlss, gamlss.dist, Rdonlp2
Additional_repositories: http://R-Forge.R-project.org
Modified: pkg/depmixS4/NAMESPACE
===================================================================
--- pkg/depmixS4/NAMESPACE 2018-09-20 12:59:46 UTC (rev 669)
+++ pkg/depmixS4/NAMESPACE 2018-10-17 09:53:34 UTC (rev 670)
@@ -10,13 +10,13 @@
importFrom(stats4, AIC, BIC, logLik, nobs, summary)
+importFrom(nlme, fdHess)
+
export(
makeDepmix,
makeMix,
lystig,
fb,
- forwardbackward,
- MVNresponse,
llratio,
multinomial,
em,
@@ -55,6 +55,7 @@
posterior,
GLMresponse,
MVNresponse,
+ forwardbackward,
transInit,
setpars,
getpars,
@@ -64,7 +65,8 @@
simulate,
summary,
logLik,
- getmodel
+ getmodel,
+ hessian
)
useDynLib(depmixS4, .registration = TRUE)
\ No newline at end of file
Modified: pkg/depmixS4/R/allGenerics.R
===================================================================
--- pkg/depmixS4/R/allGenerics.R 2018-09-20 12:59:46 UTC (rev 669)
+++ pkg/depmixS4/R/allGenerics.R 2018-10-17 09:53:34 UTC (rev 670)
@@ -54,7 +54,9 @@
setGeneric("getmodel", function(object,...) standardGeneric("getmodel"))
+setGeneric("hessian", function(object, ...) standardGeneric("hessian"))
+
# functions
setGeneric("fit", function(object, ...) standardGeneric("fit"))
Modified: pkg/depmixS4/R/fdHessHmm.R
===================================================================
--- pkg/depmixS4/R/fdHessHmm.R 2018-09-20 12:59:46 UTC (rev 669)
+++ pkg/depmixS4/R/fdHessHmm.R 2018-10-17 09:53:34 UTC (rev 670)
@@ -1,100 +1,102 @@
-
-
-hess <- function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0) {
-
- fi <- !is.null(fixed)
- cr <- !is.null(conrows)
- eq <- !is.null(equal)
-
- constr <- any(c(fi,cr,eq))
-
- if(is.nan(logLik(object))) stop("Log likelihood is 'NaN'; cannot compute hessian. ")
-
- # 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 the full set of parameters
- allpars <- getpars(object)
-
- # get the reduced set of parameters, ie the ones that will be optimized
- pars <- allpars[!fixed]
-
- 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 fit 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
- 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)
- }
- }
-
- # select only those columns of the constraint matrix that correspond to non-fixed parameters
- linconFull <- lincon
- 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]
- }
-
- # TODO: remove rows of lincon with inequality constraints!!!!
-
-
- # 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 # remove magic number here!!!!!!!!
- ans
- }
-
- return(list(fdh=fdHess(pars,logl),lin=lincon))
-}
-
-
-
-
-
-
+
+
+
+setMethod("hessian", "mix",
+ function(object,fixed=NULL,equal=NULL,
+ conrows=NULL,conrows.upper=NULL,conrows.lower=NULL,
+ method="finiteDifferences", ...) {
+
+ fi <- !is.null(fixed)
+ cr <- !is.null(conrows)
+ eq <- !is.null(equal)
+
+ constr <- any(c(fi,cr,eq))
+
+ if(is.nan(logLik(object))) stop("Log likelihood is 'NaN'; cannot compute hessian. ")
+
+ # 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 the full set of parameters
+ allpars <- getpars(object)
+
+ # get the reduced set of parameters, ie the ones that will be optimized
+ pars <- allpars[!fixed]
+
+ 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 fit 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
+ 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)
+ }
+ }
+
+ # select only those columns of the constraint matrix that correspond to non-fixed parameters
+ linconFull <- lincon
+ 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]
+ }
+
+ # TODO: remove rows of lincon with inequality constraints!!!!
+
+ # 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 # remove magic number here!!!!!!!!
+ ans
+ }
+
+ fdh <- fdHess(pars,logl)
+
+ hess <- hessian2vcov(fdh$Hessian,lincon)
+
+ return(hess)
+}
+)
\ No newline at end of file
Deleted: pkg/depmixS4/R/hess2cov.R
===================================================================
--- pkg/depmixS4/R/hess2cov.R 2018-09-20 12:59:46 UTC (rev 669)
+++ pkg/depmixS4/R/hess2cov.R 2018-10-17 09:53:34 UTC (rev 670)
@@ -1,45 +0,0 @@
-#
-# Ingmar Visser, sept 20, 2018
-#
-# Compute a corrected hessian using linear constraint matrix
-# Note: the constraints should only be the linear equality constraints, not the
-# inequality constraints!!!
-#
-
-hess2cov <- function(hess,lincon) {
- np <- dim(hess)[1]
- if(!(dim(hess)[1]==dim(hess)[2])) stop("'hess' should be a square matrix")
- se=rep(0,np)
- hs=hess
- if(nrow(lincon)>0) {
- A=lincon
- d=hs+t(A)%*%A
- di=try(solve(d),silent=TRUE)
- if(class(di)=="try-error") {
- warning("Hessian singular, ses could not be computed.")
- val=list(se=0,hs=0)
- } else {
- ada=A%*%di%*%t(A)
- adai=try(solve(ada),silent=TRUE)
- if(class(adai)=="try-error") {
- warning("Near-singular hessian, ses may be bad.\n")
- diag(ada)=diag(ada)*1.000001
- adai=try(solve(ada))
- if(class(adai)=="try-error") {
- warning("Corrected hessian also singular, ses computed without contraints.\n")
- se=sqrt(diag(di))
- } else {
- ch=di-di%*%t(A)%*%adai%*%A%*%di
- se=sqrt(diag(ch))
- }
- } else {
- ch=di-di%*%t(A)%*%adai%*%A%*%di
- se=sqrt(diag(ch))
- }
- }
- } else {
- se=sqrt(diag(solve(hs)))
- }
- val=list(se=se,hs=hs)
- val
-}
\ No newline at end of file
Copied: pkg/depmixS4/R/hessian2vcov.R (from rev 669, pkg/depmixS4/R/hess2cov.R)
===================================================================
--- pkg/depmixS4/R/hessian2vcov.R (rev 0)
+++ pkg/depmixS4/R/hessian2vcov.R 2018-10-17 09:53:34 UTC (rev 670)
@@ -0,0 +1,45 @@
+#
+# Ingmar Visser, Oct 10, 2018
+#
+# Compute a corrected covariance using linear constraint matrix and hessian
+# Note: the constraints should only be the linear equality constraints, not the
+# inequality constraints!!!
+#
+
+hessian2vcov <- function(hessian,lincon=NULL) {
+ np <- dim(hessian)[1]
+ if(!(dim(hessian)[1]==dim(hessian)[2])) stop("'hessian' should be a square matrix")
+ if(!is.null(lincon)) { # deal with linear constraints
+ nc <- ncol(lincon)
+ if(np!=nc) stop("Nr of columns in linear constraint matrix not compatible with 'hessian' dimension")
+ if(nrow(lincon)>0) {
+ A <- lincon
+ d <- hessian+t(lincon)%*%lincon
+ di <- try(solve(d),silent=TRUE)
+ if(class(di)=="try-error") {
+ warning("Hessian singular, ses could not be computed.")
+ vcov <- 0
+ } else {
+ ada <- A%*%di%*%t(A)
+ adai <- try(solve(ada),silent=TRUE)
+ if(class(adai)=="try-error") {
+ warning("Near-singular hessian, ses may be bad.\n")
+ diag(ada) <- diag(ada)*1.000001
+ adai <- try(solve(ada))
+ if(class(adai)=="try-error") {
+ warning("Corrected hessian also singular, ses computed without contraints.\n")
+ } else {
+ vcov <- di-di%*%t(A)%*%adai%*%A%*%di
+ }
+ } else {
+ vcov <- di-di%*%t(A)%*%adai%*%A%*%di
+ }
+ }
+ } else { # linear constraint matrix present but 0 rows
+ vcov <- solve(hessian)
+ }
+ } else { # no linear constraints
+ vcov <- solve(hessian)
+ }
+ vcov
+}
More information about the depmix-commits
mailing list