[Depmix-commits] r673 - pkg/depmixS4/R

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

Author: ingmarvisser
Date: 2018-10-17 17:50:42 +0200 (Wed, 17 Oct 2018)
New Revision: 673

Modified:
   pkg/depmixS4/R/hessian.R
Log:
=finished hessian function, now includes checks for boundary parameters

Modified: pkg/depmixS4/R/hessian.R
===================================================================
--- pkg/depmixS4/R/hessian.R	2018-10-17 14:12:03 UTC (rev 672)
+++ pkg/depmixS4/R/hessian.R	2018-10-17 15:50:42 UTC (rev 673)
@@ -79,7 +79,7 @@
 		lin.l <- c(lin.l,rep(0,nrow(equal)))				
 	}
 	
-	# incorporate general linear constraints, if any
+	# 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)
@@ -100,26 +100,24 @@
 	# 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)
 	
-	print(low)
-	print(up)
 	
-	print(which(fixed))
+	if(length(which(fixed)>0)) elements[which(fixed)] <- "fix"
+	if(length(bnd)>0) elements[bnd] <- "bnd"
 	
-	print(union(low, up, fixed))
-	
 	# get the reduced set of parameters, ie the ones that the hessian will be computed for
 	# only non-fixed parameters
-	pars <- allpars[!fixed]
+	pars <- allpars[which(elements=="inc")]
 	
 	# select only those columns of the constraint matrix that correspond to non-fixed parameters
-	linconFull <- lincon
-	lincon <- lincon[,!fixed,drop=FALSE]
+	lincon <- lincon[,which(elements=="inc"),drop=FALSE]
 	
 	# remove redundant rows in lincon (all zeroes)
 	allzero <- which(apply(lincon,1,function(y) all(y==0)))
@@ -129,15 +127,21 @@
 		lin.l <- lin.l[-allzero]
 	}
 	
-	# TODO: remove rows of lincon with inequality constraints!!!!
+	# 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]
+	}
 	
-	
 	# make loglike function that only depends on pars
 	logl <- function(pars) {
-		allpars[!fixed] <- pars
+		allpars[which(elements=="inc")] <- pars
 		object <- setpars(object,allpars)
 		ans = -as.numeric(logLik(object))
-		if(is.na(ans)) ans = 100000 # remove magic number here!!!!!!!!
+		if(is.na(ans)) ans = 1000000 # remove magic number here!!!
 		ans
 	}