Permalink
Cannot retrieve contributors at this time
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
74 lines (73 sloc)
2.98 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
\*********************************************************************************/ | |
\* This is the Gill/Murray cholesky routine. 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 R format. */ | |
\*********************************************************************************/ | |
proc gmchol(A); | |
/* calculates the Gill-Murray generalized choleski decomposition */ | |
/* input matrix A must be non-singular and symmetric */ | |
/* Author: Jeff Gill. Part of the Hessian Project. */ | |
local i,j,k,n,sum,R,theta_j,norm_A,phi_j,delta,xi_j,gamm,E,beta_j; | |
n = rows(A); | |
R = eye(n); | |
E = zeros(n,n); | |
norm_A = maxc(sumc(abs(A))); | |
gamm = maxc(abs(diag(A))); | |
delta = maxc(maxc(__macheps*norm_A~__macheps)); | |
for j (1, n, 1); | |
theta_j = 0; | |
for i (1, n, 1); | |
sum = 0; | |
for k (1, (i-1), 1); | |
sum = sum + R[k,i]*R[k,j]; | |
endfor; | |
R[i,j] = (A[i,j] - sum)/R[i,i]; | |
if (A[i,j] -sum) > theta_j; | |
theta_j = A[i,j] - sum; | |
endif; | |
if i > j; | |
R[i,j] = 0; | |
endif; | |
endfor; | |
sum = 0; | |
for k (1, (j-1), 1); | |
sum = sum + R[k,j]^2; | |
endfor; | |
phi_j = A[j,j] - sum; | |
if (j+1) <= n; | |
xi_j = maxc(abs(A[(j+1):n,j])); | |
else; | |
xi_j = maxc(abs(A[n,j])); | |
endif; | |
beta_j = sqrt(maxc(maxc(gamm~(xi_j/n)~__macheps))); | |
if delta >= maxc(abs(phi_j)~((theta_j^2)/(beta_j^2))); | |
E[j,j] = delta - phi_j; | |
elseif abs(phi_j) >= maxc(((delta^2)/(beta_j^2))~delta); | |
E[j,j] = abs(phi_j) - phi_j; | |
elseif ((theta_j^2)/(beta_j^2)) >= maxc(delta~abs(phi_j)); | |
E[j,j] = ((theta_j^2)/(beta_j^2)) - phi_j; | |
endif; | |
R[j,j] = sqrt(A[j,j] - sum + E[j,j]); | |
endfor; | |
retp(R'R); | |
endp; |