[Depmix-commits] r641 - in pkg/depmixS4: R man

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

Author: ingmarvisser
Date: 2014-10-09 16:44:58 +0200 (Thu, 09 Oct 2014)
New Revision: 641

Added:
   pkg/depmixS4/R/hess2cov.R
Modified:
   pkg/depmixS4/man/depmix-internal.Rd
Log:
Added function hess2cov that computes the corrected hessian given an (approximate) hessian subject to linear equality constraints; the constraint matrix is passed as argument next the hessian.

Added: pkg/depmixS4/R/hess2cov.R
===================================================================
--- pkg/depmixS4/R/hess2cov.R	                        (rev 0)
+++ pkg/depmixS4/R/hess2cov.R	2014-10-09 14:44:58 UTC (rev 641)
@@ -0,0 +1,42 @@
+# 
+# Ingmar Visser, oct 9, 2014
+# 
+# Compute a corrected hessian using
+# 
+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


Property changes on: pkg/depmixS4/R/hess2cov.R
___________________________________________________________________
Added: svn:eol-style
   + native

Modified: pkg/depmixS4/man/depmix-internal.Rd
===================================================================
--- pkg/depmixS4/man/depmix-internal.Rd	2014-09-29 19:05:02 UTC (rev 640)
+++ pkg/depmixS4/man/depmix-internal.Rd	2014-10-09 14:44:58 UTC (rev 641)
@@ -10,6 +10,8 @@
 \alias{makeResponseModels}
 \alias{makeTransModels}
 
+\alias{hess2cov}
+
 \alias{mlogit}
 \alias{multinomial}