; ======================================================================
;
; Function ConfInt
;
; Calculate uncertainty on model parameters which were fit by
; minimization of chisquare.
;
; Written By:  BA Franz, ARC, 2/95
;
; Calling Sequence:
;     sigma = ConfInt(covar,nfree,chisqr,confidence=confidence) 
;
; Inputs:
;     covar  - m x m covariance matrix (or hessian) of fitted parameters 
;     nfree  - number of degrees of freedom (e.g., n - m)
;     chisqr - reduced chisquare of the fit (i.e., per degree of freedom)
;
; Optional Keyword Inputs:
;     confidence - desired confidence level (def=68.0)
;     hessian    - set this to 1 if hessian was supplied in place 
;                  the covariance matrix
;
; Output:
;     sigma - m-element vector of parameter uncertainties
;
; Description of Algorithm (by T. J. Sodroski)
; The parameter uncertainties are obtained from an m-dimensional
; ellipsoidal confidence region in m-parameter space.  The length
; of each principal axis i of the confidence region is given by
; sqrt(2E/Li), where E is the largest allowable deviation of
; chisquare from its optimal value for a given level of significance
; (as derived from the F-distribution), and Li is the associated 
; eigenvalue of the hessian matrix for parameter i.  The components
; of each principal axis, Pi, in parameter space are then obtained
; by: [Pi] = [U] * [Wi], where [U] is the unitary transformation
; matrix whose columns are the normalized eigenvectors of the
; Hessian, and [Wi] is the principal axis vector in the coordinate
; system defined by the axes of the ellipsoid.  The uncertainty in
; each parameter is then given by the maximum component of the 
; principal axis vectors along the parameter axis.
;
; References: 
; Bard, Y.  1974, Nonlinear Parameter Estimation, (Orlando: Academic)
; Sodroski, T. J.  1988, Ph.D. thesis, University of Maryland
;
; ======================================================================
function ConfInt,covar,nfree,chisqr,confidence=conf,hessian=hessian 

if (n_elements(conf) eq 0) then conf=68.0

;
; Get number of parameters from size of COVAR
;
s = size(covar)
if ( (s(0) ne 2) or (s(1) ne s(2)) ) then begin
    message,'Covariance matrix must be a square array',/info
   return,-1
endif else nterms = s(1)
diag = indgen(nterms)*(nterms+1)         ; Subscripts of diagonal elements

;
; Get Hessian from COVAR unless Hessian was passed instead
;
if (not keyword_set(Hessian)) then begin
    hess = invert(covar,status)
    case (status) of
      0:
      1: message,'Inversion of covariance matrix failed',/info
      2: message,'Inversion of covariance matrix inaccurate',/info
    endcase
endif else hess = covar

;
; Use f-distribution to determine ep, the largest allowable deviation
; of chisquare from its optimal value within the confidence region.
;
fd = f_test(1.-conf/100.,nterms,nfree)
ep = 2.0 * nterms * chisqr * fd

;
; Compute confidence intervals for each paramter versus other parameters
;
temp1  = hess                             ; The unnormalized hessian
nr_tred2,temp1,temp2,temp3,/double        ; Reduce to tridiagonal form
nr_tqli,temp2,temp3,temp1,/double         ; Compute:
eig    = temp2                            ;   Eigenevalues
eigvec = temp1                            ;   Eigenvectors
paxis_len = dblarr(nterms,nterms)         ; Length of principal axes of
paxis_len(diag) = sqrt((ep/eig) > 1e-30)  ;   confidence region
theta  = abs(paxis_len#transpose(eigvec)) ; Projection of principal axes of 
                                          ;   conf regions onto param axes
;
; Compute parameter uncertainties as max components of principal axes
; in parameter space.
;
sigma = fltarr(nterms)                    
for i=0,nterms-1 do sigma(i) = max(theta(*,i))            

return,sigma
end