[Depmix-commits] r669 - in pkg/depmixS4: . R

[noreply at r-forge.r-project.org]

Author: ingmarvisser
Date: 2018-09-20 14:59:46 +0200 (Thu, 20 Sep 2018)
New Revision: 669

Added:
   pkg/depmixS4/R/fdHessHmm.R
Modified:
   pkg/depmixS4/DESCRIPTION
   pkg/depmixS4/NEWS
   pkg/depmixS4/R/hess2cov.R
Log:
=added fdHessHmm.R

Modified: pkg/depmixS4/DESCRIPTION
===================================================================
--- pkg/depmixS4/DESCRIPTION	2018-09-11 15:03:34 UTC (rev 668)
+++ pkg/depmixS4/DESCRIPTION	2018-09-20 12:59:46 UTC (rev 669)
@@ -1,6 +1,6 @@
 Package: depmixS4
-Version: 1.3-5
-Date: 2018-08-24 (R666)
+Version: 1.4-0
+Date: 2018-09-20
 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>

Modified: pkg/depmixS4/NEWS
===================================================================
--- pkg/depmixS4/NEWS	2018-09-11 15:03:34 UTC (rev 668)
+++ pkg/depmixS4/NEWS	2018-09-20 12:59:46 UTC (rev 669)
@@ -1,3 +1,15 @@
+Changes in depmixS4 version 1.4-0
+
+Major changes
+
+  o Added functionality for computing a finite differences hessian of fitted 
+    models, which is returned by default from the fit function. 
+	
+Minor changes
+
+  o Multivariate normal distribution models now use ML estimates (instead of the 
+    unbiased version that was hitherto used). 	
+
 Changes in depmixS4 version 1.3-5
 
   o GLM responses now allow missing values in covariates, as long as the

Added: pkg/depmixS4/R/fdHessHmm.R
===================================================================
--- pkg/depmixS4/R/fdHessHmm.R	                        (rev 0)
+++ pkg/depmixS4/R/fdHessHmm.R	2018-09-20 12:59:46 UTC (rev 669)
@@ -0,0 +1,100 @@
+
+
+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))
+}
+
+
+
+
+
+


Property changes on: pkg/depmixS4/R/fdHessHmm.R
___________________________________________________________________
Added: svn:executable
   + *

Modified: pkg/depmixS4/R/hess2cov.R
===================================================================
--- pkg/depmixS4/R/hess2cov.R	2018-09-11 15:03:34 UTC (rev 668)
+++ pkg/depmixS4/R/hess2cov.R	2018-09-20 12:59:46 UTC (rev 669)
@@ -1,8 +1,11 @@
 # 
-# Ingmar Visser, oct 9, 2014
+# Ingmar Visser, sept 20, 2018
 # 
-# Compute a corrected hessian using
+# 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")