From 71b1b000b3c8b20dd7bd5b85c97346c560da2457 Mon Sep 17 00:00:00 2001 From: Dan H Date: Wed, 14 Jun 2017 22:45:56 -0400 Subject: [PATCH] Add important fix to parameter error. Numeric hessians from numdifftools are not well behaved at the moment (for both chi^2 and the fit function). Tolerances either need to be increased, the second derivative(s) of the fit function need(s) to be set to 0, or the Hessian for the fit function should be computed analytically. TODO: Add an option to enter the first derivatives of the fit function by hand, GEVP analysis, and double jackknife resampling. This commit should be nearly at the 1.0 release stage, and should be considered production ready. --- latfit/finalout/geterr.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/latfit/finalout/geterr.py b/latfit/finalout/geterr.py index 13ebc21c..a288f260 100644 --- a/latfit/finalout/geterr.py +++ b/latfit/finalout/geterr.py @@ -14,29 +14,30 @@ def geterr(result_min, covinv, coords): param_err = ["err" for i in range(len(result_min.x))] if result_min.fun > 0 and result_min.status == 0: - HFUNC = lambda xrray: chi_sq(xrray, covinv, coords) - HFUN = nd.Hessian(HFUNC) + #Numeric Hessian of chi^2 is unstable. Use at own risk. + #HFUNC = lambda xrray: chi_sq(xrray, covinv, coords) + #HFUN = nd.Hessian(HFUNC) flag = 0 cutoff = 1/(CUTOFF*SCALE) l_coords=len(coords) l_params=len(result_min.x) - #delete this! - if len(HFUN(result_min.x)) != len(result_min.x): - print("huh.") - sys.exit(1) - #compute the gradient and hessian of the fit function def g(x): return np.array([fit_func(coords[i][0],x) for i in range(len(coords))]) - fhess=nd.Gradient(nd.Gradient(g))(result_min.x) + #fhess=nd.Gradient(nd.Gradient(g))(result_min.x) + #TODO: allow for this to be entered by hand. grad=nd.Gradient(g)(result_min.x) #compute analytically hessian of chi^2 with respect to fit parameters if l_params!=1: - Hess=[[2*fsum(np.dot(fhess[i][a][b]*covinv[i][j],(g(result_min.x)[j]-coords[j][1]))+np.dot(grad[i][a]*covinv[i][j],grad[j][b]) for i in range(l_coords) for j in range(l_coords)) for b in range(l_params)] for a in range(l_params)] + #use numeric second deriv of fit fit func. also unstable (and often wrong) + #Hess=[[2*fsum(np.dot(fhess[i][a][b]*covinv[i][j],(g(result_min.x)[j]-coords[j][1]))+np.dot(grad[i][a]*covinv[i][j],grad[j][b]) for i in range(l_coords) for j in range(l_coords)) for b in range(l_params)] for a in range(l_params)] + Hess=[[2*fsum(np.dot(grad[i][a]*covinv[i][j],grad[j][b]) for i in range(l_coords) for j in range(l_coords)) for b in range(l_params)] for a in range(l_params)] else: - Hess=[[2*fsum(np.dot(fhess[i]*covinv[i][j],(g(result_min.x)[j]-coords[j][1]))+np.dot(grad[i]*covinv[i][j],grad[j]) for i in range(l_coords) for j in range(l_coords))]] + #see above about instability of fit func hessian. this else block is for one parameter fits. + #Hess=[[2*fsum(np.dot(fhess[i]*covinv[i][j],(g(result_min.x)[j]-coords[j][1]))+np.dot(grad[i]*covinv[i][j],grad[j]) for i in range(l_coords) for j in range(l_coords))]] + Hess=[[2*fsum(np.dot(grad[i]*covinv[i][j],grad[j]) for i in range(l_coords) for j in range(l_coords))]] #compute hessian inverse of chi^2 try: @@ -67,6 +68,7 @@ def g(x): print("leading to complex errors in some or all of the") print("fit parameters. failing entry:",2*HINV[i][i]) print("Attempting to continue with entries zero'd below", cutoff) + print("The errors in fit parameters which result shouldn't be taken too seriously.") print('Debugging info:') print('2x Hessian Inverse:') print(2*HINV)