forum.alglib.nethttp://forum.alglib.net/ Accessing Hessian or unscaled covariance of LM fit?http://forum.alglib.net/viewtopic.php?f=2&t=4601 Page 1 of 1

 Author: elivesay2000 [ Fri Apr 05, 2024 8:40 pm ] Post subject: Accessing Hessian or unscaled covariance of LM fit? Hello, I have a question regarding getting access to the Hessian or the unscaled covariance following a Levenberg Marqurdt fit. I have the functional form of the equation as well as the derivatives of that functional form with respect to the model parameters. However, I do not have an analytial expression for the Hessian. As a result, I am using lsfitcreatefg (as opposed to the related "f" or fgh" versions of lsfitcreate) like this: double epsx = 0.000001; int maxits = 5000; int info = 0; alglib.lsfitstate state; alglib.lsfitreport report; alglib.lsfitcreatefg (x, y, c, false, out state); alglib.lsfitsetcond (state, epsx, maxits); alglib.lsfitfit (state, function_myEq_func, function_myEq_grad, null, null); alglib.lsfitresults (state, out info, out c, out report); The fit works really well. It returns the same coefficients that I get when I perform the fit with Python. That said, how can I obtain the Hessian evaluated at the minimizing point? It looks like the Hessian should be available as state.innerobject.h. Although h is the nomenclature for the Hessian in the calls to lsfit (as noted above) and h is a 2D array of doubles, h has a size of 0x0. Is there a flag that I need to set to ask that the Hessian be reported? Related to this is the covariance matrix that is reported in report.covpar. This value is equivalent to the unscaled covariance matrix that is reported by a call to the Python Method minimize with the parameter scale_covar = true (its default value). With scale_covar = false, the unscaled covariance matrix is returned by Python. How can I either "unscale" the covariance matrix reported as report.covpar or ask ALGLIB to return the unscaled version of the covariance matrix?This question looks similar to this one: http://forum.alglib.net/viewtopic.php?f=2&t=293&hilit=covariance

 Author: elivesay2000 [ Mon Apr 08, 2024 2:58 am ] Post subject: Re: Accessing Hessian or unscaled covariance of LM fit? I'm following-up here with an implemented example.In the image below is an image of functioning example of Python. For a single dataset, a fit to a second order polynomial (y = a0 + a1*x + a2*x^2) is performed and either the scaled or unscaled covariance matrix is returned.Attachment: Python_scale_covar_question.png [ 122.29 KiB | Viewed 2492 times ] Below is code that utilizes ALGLIB to perform a fit to the same functional form (i.e., a second order polynomial) of the same dataset. The same coefficients are found. Based on a check of the diagonal elements (see the output at the bottom of the image), the covariance matrix that is calculated by ALGLIB is the scaled covariance matrix - the diagonal values agree with the diagonal elements of the scaled covariance matrix returned by Python, not the diagonal elements of the unscaled covariance matrix returned by Python.How can I get the unscaled covariance matrix using ALGLIB?Additionally, how can I get the Hessian evaluated at the minimizing point?Code:      [TestMethod]      public void TestAlgLibSecondOrderModel()      {         var x      = new double[,] { { 0 },          { 1 },               { 2 },           { 3 },            { 4 }             };        //var y      = new double[]  { 3.9 + 0.0 + 0.0, 3.9 + 2.0 + (-0.0), 3.9 + 8.0 + 0.0, 3.9 + 18.0 + 0.0, 3.9 + 32.0 + -0.0 };   // No noise          var y      = new double[]  { 3.9 + 0.0 + 0.3, 3.9 + 2.0 + (-0.2), 3.9 + 8.0 + 0.4, 3.9 + 18.0 + 0.5, 3.9 + 32.0 + -0.4 };         var c      = new double[]  { 10.0, 0.0, 0.0 };         var epsx   = 0.000001;         var maxits   = 5000;         var info   = 0;         alglib.lsfitstate state;         alglib.lsfitreport rep;         // Fitting without weights         alglib.lsfitcreatefg   (x, y, c, false, out state);         alglib.lsfitsetcond      (state, epsx, maxits);         alglib.lsfitfit         (state, function_SecondOrder_func, function_SecondOrder_grad, null, null);         alglib.lsfitresults      (state, out info, out c, out rep);         // Assertions         // Expected values come from Python & are also found in TestMathnetLMSolver.TestLMSolverResultsForEquation25aFitUnconstrainedCoefficients()         Assert.AreEqual(3.97, c[0], 0.01);         Assert.AreEqual(0.30, c[1], 0.01);         Assert.AreEqual(1.91, c[2], 0.01);                 // Check the solver's scaled covariance results         var cov = rep.covpar;         Assert.AreEqual(0.20295512, cov[0, 0], 1.00);  // From Python (these are scaled covariance values).         Assert.AreEqual(0.28479186, cov[1, 1], 0.01);  // From Python (these are scaled covariance values).         Assert.AreEqual(0.01636735, cov[2, 2], 0.01);  // From Python (these are scaled covariance values).      }For completeness, here is the Python code:Code:from lmfit   import minimize, Parametersimport numpy as np# The Method that encapsulates the second order polynomial Model that will be fit to measured pressure datadef SecondOrderModel(params, x):    a0      = params['a0']    a1      = params['a1']    a2      = params['a2']    returnValue = a0 + a1 * x + a2 * pow(x, 2)    return(returnValue)# This Method returns the difference between the measured and fitted values and is needed by the LM fitting routinedef residualSecondOrderModel(params, x, y, uncertainty):    theModel    = SecondOrderModel(params, x)    returnValue = (y - theModel)/(uncertainty)    return(returnValue)# --- Main Routine ---# Initial Guesses (IGs).IG_a0 = 10.0IG_a1 =  0.0IG_a2 =  0.0x_ToFit = np.array([0.0,             1.0,                2.0,             3.0,              4.0              ])#y_ToFit = np.array([3.9 + 0.0 + 0.0, 3.9 + 2.0 + (-0.0), 3.9 + 8.0 + 0.0, 3.9 + 18.0 + 0.0, 3.9 + 32.0 + -0.0]) # 3.9 + 0x + 2x^2 with no noisey_ToFit = np.array([3.9 + 0.0 + 0.3, 3.9 + 2.0 + (-0.2), 3.9 + 8.0 + 0.4, 3.9 + 18.0 + 0.5, 3.9 + 32.0 + -0.4]) # 3.9 + 0x + 2x^2 with a bit of hand-made noise# Set all uncertainties to 1.0 - Numerical Recipes, Third Edition, Section 15.4 (page 788) & initialize parametersuncertainty     = np.ones(x_ToFit.size)params2ndOrder  = Parameters()params2ndOrder.add("a0",       value = IG_a0)params2ndOrder.add("a1",       value = IG_a1)params2ndOrder.add("a2",       value = IG_a2)out2ndOrder_lmfit_T = minimize(residualSecondOrderModel, params2ndOrder, args = (x_ToFit, y_ToFit, uncertainty), scale_covar=True)print("--- SCALED COVARIANCE (scale_covar = True) ---")print(out2ndOrder_lmfit_T.params["a0"])print(out2ndOrder_lmfit_T.params["a1"])print(out2ndOrder_lmfit_T.params["a2"])print(out2ndOrder_lmfit_T.covar)out2ndOrder_lmfit_F = minimize(residualSecondOrderModel, params2ndOrder, args = (x_ToFit, y_ToFit, uncertainty), scale_covar=False)print("--- UNSCALED COVARIANCE (scale_covar = False) ---")print(str(out2ndOrder_lmfit_F.params["a0"]) + "   a0_T/a0_F = " + str(out2ndOrder_lmfit_T.covar[0][0] / out2ndOrder_lmfit_F.covar[0][0]))print(str(out2ndOrder_lmfit_F.params["a1"]) + "   a1_T/a1_F = " + str(out2ndOrder_lmfit_T.covar[1][1] / out2ndOrder_lmfit_F.covar[1][1]))print(str(out2ndOrder_lmfit_F.params["a2"]) + "   a2_T/a2_F = " + str(out2ndOrder_lmfit_T.covar[2][2] / out2ndOrder_lmfit_F.covar[2][2]))print(out2ndOrder_lmfit_F.covar)

 Author: elivesay2000 [ Mon Apr 08, 2024 9:16 pm ] Post subject: Re: Accessing Hessian or unscaled covariance of LM fit? After a bit more investigation I've found that the scaling factor is the reduced chi squared statistic. The rep.noise variable is an array of doubles whose size is dictated by the number of elements in the dataset. In my example, I have 5 data points (one each x = 0, 1, 2, 3, 4) and the square root of the chi squared statistic is found in the rep.noise variable, all 5 elements of which are equal to 0.47868868499564171. The square of this value is 0.229142857 which is the value of the chi squared statistic reported by Python.This begs the question of why is the chi squared statistic reported once per data point?My C# code that utilizes ALGLIB has been updated to reflect this scaling factor (see below).Code:      [TestMethod]      public void TestAlgLibSecondOrderModel()      {         var x      = new double[,] { { 0 },          { 1 },               { 2 },           { 3 },            { 4 }             };       //var y      = new double[]  { 3.9 + 0.0 + 0.0, 3.9 + 2.0 + (-0.0), 3.9 + 8.0 + 0.0, 3.9 + 18.0 + 0.0, 3.9 + 32.0 + -0.0 };   // No noise         var y      = new double[]  { 3.9 + 0.0 + 0.3, 3.9 + 2.0 + (-0.2), 3.9 + 8.0 + 0.4, 3.9 + 18.0 + 0.5, 3.9 + 32.0 + -0.4 };         var c      = new double[]  { 10.0, 0.0, 0.0 };         var epsx   = 0.000001;         var maxits   = 5000;         var info   = 0;         alglib.lsfitstate state;         alglib.lsfitreport rep;         // Fitting without weights         alglib.lsfitcreatefg   (x, y, c, false, out state);         alglib.lsfitsetcond      (state, epsx, maxits);         alglib.lsfitfit         (state, function_SecondOrder_func, function_SecondOrder_grad, null, null);         alglib.lsfitresults      (state, out info, out c, out rep);         // Assertions         // Expected values come from Python & are also found in TestMathnetLMSolver.TestLMSolverResultsForEquation25aFitUnconstrainedCoefficients()         Assert.AreEqual(3.97, c[0], 0.01);         Assert.AreEqual(0.30, c[1], 0.01);         Assert.AreEqual(1.91, c[2], 0.01);                 // Check the solver's scaled covariance results         var cov = rep.covpar;         Assert.AreEqual(0.20295512, cov[0, 0], 0.001);  // From Python (these are scaled covariance values).         Assert.AreEqual(0.28479186, cov[1, 1], 0.001);  // From Python (these are scaled covariance values).         Assert.AreEqual(0.01636735, cov[2, 2], 0.001);  // From Python (these are scaled covariance values).         // Check the solver's unscaled covariance results         var scalingFactor = Math.Pow(rep.noise[0], 2);   // Noise has 5 elements (the same number of data points defined above) that all appear to be identical. Will the noise in every cell always be identical?         Assert.AreEqual(0.88571438, cov[0, 0] / scalingFactor, 0.001);  // From Python (these are unscaled covariance values).         Assert.AreEqual(1.24285726, cov[1, 1] / scalingFactor, 0.001);  // From Python (these are unscaled covariance values).         Assert.AreEqual(0.07142857, cov[2, 2] / scalingFactor, 0.001);  // From Python (these are unscaled covariance values).      }Aside - the following Python article on github is a useful reference for the scaled covariance matrix.Attachment: 2024-04-08 13_30_39-lmfit-py_doc_fitting_rst at master_lmfit_lmfit-py.png [ 153.05 KiB | Viewed 2489 times ]

 Author: Sergey.Bochkanov [ Sat Apr 13, 2024 3:11 pm ] Post subject: Re: Accessing Hessian or unscaled covariance of LM fit? Hi!Sorry for the delay in replying! The lsfitstate.h field is intended for communication with the user and this particular protocol is no longer publicly supported, so you can't get Hessian from it. It is possible to include Hessian into reports, I will add it to the issues tracker. However, I am not sure whether it will be 4.02 which is scheduled for April-May, or 4.03. I will reply with my decision later.

 Author: elivesay2000 [ Mon Apr 15, 2024 3:46 am ] Post subject: Re: Accessing Hessian or unscaled covariance of LM fit? Understood. Thanks for the response!Related to this I have been looking at the noise field of the report class. In the simple example that I described in my April 08 post I was fitting 5 data points and the noise field was an array of doubles that had 5 elements and all 5 elements were equal to the same value, namely 0.47868868499564171. Is it fair to assume that if I fit N data points the noise field will be an array of N doubles and all N of them will be equal to the same value?

 Page 1 of 1 All times are UTC Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Grouphttp://www.phpbb.com/