################################################################################# # This is a set of test routines for the sechol.s function. Reference: # # # # Gill, Jeff and Gary King. ``What to do When Your Hessian is Not Invertible: # # Alternatives to Model Respecification in Nonlinear Estimation,'' Sociological # # Methods and Research, Vol. 32, No. 1 (2004): Pp. 54--87. # # # # Abstract # # # # What should a researcher do when statistical analysis software terminates # # before completion with a message that the Hessian is not invertable? The # # standard textbook advice is to respecify the model, but this is another way # # of saying that the researcher should change the question being asked. # # Obviously, however, computer programs should not be in the business of # # deciding what questions are worthy of study. Although noninvertable # # Hessians are sometimes signals of poorly posed questions, nonsensical # # models, or inappropriate estimators, they also frequently occur when # # information about the quantities of interest exists in the data, through # # the likelihood function. We explain the problem in some detail and lay out # # two preliminary proposals for ways of dealing with noninvertable Hessians # # without changing the question asked. # # # # Also available is the software to implement the procedure described in this # # paper in Gauss format. # ################################################################################# # FIRST EXAMPLE X <- matrix(c(5,2,3,2.95,2,1,2,1,5,2,3,3),ncol=3) Y <- c(9,11,-5,-2) b.hat <- solve(t(X)%*%X)%*%t(X)%*%Y Y.hat <- X%*%b.hat e <- Y - Y.hat sigma <- as.numeric(t(e)%*%e) b.var <- sigma/(nrow(X)-ncol(X)) * solve(t(X)%*%X) b.cor <- b.var for (i in 1:nrow(b.cor)) { for (j in 1:ncol(b.cor)) { if (i != j) b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j]) else b.cor[i,j] <- sqrt(b.var[i,j]) } } # SECOND EXAMPLE X <- matrix(c(5,2,3,2.99,2,1,2,1,5,2,3,3),ncol=3) Y <- c(9,11,-5,-2) b.hat <- solve(t(X)%*%X)%*%t(X)%*%Y Y.hat <- X%*%b.hat e <- Y - Y.hat sigma <- as.numeric(t(e)%*%e) b.var <- sigma/(nrow(X)-ncol(X)) * solve(t(X)%*%X) b.cor <- b.var for (i in 1:nrow(b.cor)) { for (j in 1:ncol(b.cor)) { if (i != j) b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j]) else b.cor[i,j] <- sqrt(b.var[i,j]) } } ginv <- function(X, tol=sqrt(.Machine$double.eps)) { s <- svd(X) nz <- s$d > tol * s$d[1] if(any(nz)) s$v[,nz] %*% (t(s$u[, nz])/s$d[nz]) else X*0 } mpinv <- function(X, tol = sqrt(.Machine$double.eps)) { s <- svd(X) e <- s$d e[e > tol] <- 1/e[e > tol] s$v %*% diag(e) %*% t(s$u)) } # THIRD EXAMPLE X <- matrix(c(5,2,3,2.999,2,1,2,1,5,2,3,3),ncol=3) Y <- c(9,11,-5,-2) b.hat <- mpinv(t(X)%*%X)%*%t(X)%*%Y Y.hat <- X%*%b.hat e <- Y - Y.hat sigma <- as.numeric(t(e)%*%e) b.var <- sigma * mpinv(t(X)%*%X) b.cor <- b.var for (i in 1:nrow(b.cor)) { for (j in 1:ncol(b.cor)) b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j]) } diag(b.cor) <- sqrt(diag(b.var)) # FIRST CHOLESKY EXAMPLE S <- matrix(c(2,0,2.4,0,2,0,2.4,0,3),ncol=3) chol(S) T <- sechol(S) t(T) %*% T # SECOND CHOLESKY EXAMPLE S <- matrix(c(2,0,2.5,0,2,0,2.5,0,3),ncol=3) sechol(S) # THIRD CHOLESKY EXAMPLE S <- matrix(c(2,0,10,0,2,0,10,0,3),ncol=3) T <- sechol(S) t(T) %*% T