Astro-Stats & Python : Levenberg-Marquardt Statistics

Another exciting day in my first course in Grad level Astro-Statistics.  Today I'll discuss a simple exercise that was great in helping me comprehend and apply Levenberg-Marquardt Statistics.

As always, you can follow along and run my code if such is your heart's desire.  However, this time I wrote code in an IPython Notebook.  You can just look at it on my Github page here, or if you have Anaconda installed onto your computer you can play with it yourself if you copy the file from here and paste into your own script.

As usual, I must first simulate data as if I collected them from some star or astrophysical source.  I did this with a random number generator (of a Gaussian distribution) in python which acted as my "noise".  From an astronomer's perspective, noise is just unwanted, additional signal/light from somewhere in space or on Earth that was measured by the telescope's camera.  When observing a star in outer space, you only want your telescope to receive the light from that specific star.  This does not happen because of how the earth's atmosphere alters and adds on to the measurements of that star's light as that light travels from outer space to earth's atmosphere to the Earth's surface.

Figure 1: Plot of simulated data.  Just by looking at the plot, you can see why the shape of a plot that follows a Gaussian distribution is commonly referred as a bell curve.

 In previous posts, I generated data by using a polynomial equation where the user determined how many coefficients were involved and thus how many terms where in that equation.  Subsequently, a random number would be added to the equation to make it more realistic and unpredictable.  For this exercise I use a different equation that specifically follows a Gaussian distribution.  Furthermore, at the end of this equation I still added noise by using a random number generator.  So now, I have data that abides by a Gaussian along with the confidence interval (or uncertainties) of the measured data points abiding by a Gaussian as well.  My Gaussian equation is

\[G(x_i) = A e^{-\frac{1}{2}\Big( \frac{B-x_i}{C} \Big)^2 } + E  \]

where G is a function of the x-axis values \(x_i\), A is for the amplitude, B is the centroid of the Gaussian, C is the width and E is for the offset.  A  and E are known as free parameters and are usually insignificant when analyzing astrophysical data.  This equation will also serve as my model for the Levenberg-Marquardt (LM) component of this exercise.

Now that I have my simulated data and I've chosen my model, we must understand how to perform the least-squares fitting on this data.  In this case, I must find the least \(\chi^2\) in order to find the best fitting parameters A, B, C and E for my data.  The equation is

\[ \chi^2 = \sum_{i}^n \Big(\frac{D_i - G_i}{\sigma_{D_i}}   \Big)^2  \]

where \(D_i\) is the measurement (or data point), \(G_i\) is the model's result and \(\sigma\) is the uncertainty of the measurement \(D_i\).  There is a standard python function that I use to find which combination of values for the 4 parameters will yield a \(G_i\) that produces the least \(\chi^2\).  That function is called "leastsq." However, I first must create a function that would input the residual into "leastsq" function:

\[ R = \frac{D_i - G_i}{\sigma_{D_i}} \]
where R  is the residual.

In order to begin utilizing these equations, I must make initial guesses as to what the value of the parameters are.  This will result in the first \(G_i\) calculated and thus the first residual and \(\chi^2\) will be found.  From there, the program will continue to try different values and gravitate towards values that yield a smaller and smaller \(\chi^2\).

For this example, I guessed that the parameters were A = 20, B = 7, C = 5 and E = 9.   The true values of the parameters were A = 10, B = 5, C = 2 and E = 4.

After finding the best values for my 4 fitting parameters, I plug those values into my model and check to see how well this specific model aligns with the data.

Figure 2: Top- my code showing the values that yielded the least \(\chi^2\), which consequently fits the data better than any other combination of values tried.  Bottom- Plot of original data and the fit for that data.  The fit seems to be incredibly agreeable with the data.

In summary, in order to use the powerful statistical method derived by Levenberg and Marquardt, I had to make initial guesses for my parameters, use the \(\chi^2\) equation and find the combination of parameter-values that yielded the least \(\chi^2\).  LM statistics allowed me to find the best fitting parameters for my data.  But it did not provide the confidence interval of those parameters.   Knowing the uncertainty of the evaluated parameters is very important!  This is where Monte Carlo (MC) and Bootstrapping statistics can save the day.  If you ever have a complex and large set of data for which you must find a satisfactory fit to, I recommend that you use LM statistics to find the best fitting parameter first and then use MC and Bootstrapping statistics to find the confidence interval of each parameter.