No uncertainties being returned with LevMar fitter

[karllark]

As noted in PRs #141, #142, #143, we are not getting uncertainties when using the LevMar fitter interface provided by in astropy.modeling. Initial thoughts were that the problem could be the Compound astropy.model we use or that some of the components in this Compound model were not "stock" astropy functions.

Work by @alexmaragko showed that it was not the Compound model.

Continuing this investigation, I have been able to show that it is not the use of non-"stock" models. The issue seems to be when a parameter is not needed or is at one of the limits (zero in the example 2 below). This is not an issue with astropy.modeling, but is the result of the underlying scipy function not returning a covariance matrix. Astropy.modeling uses the scipy returned covariance matrix to calculate the uncertainties.

1st example works - uncertainties reported. 2nd example does not report uncertainties. Only difference is that the simulation has tau_sil = 0 instead of tau_sil = 1.

Example 1 of fitting a Compound model with (Drude + line) * S07_attenuation model (see code below). The tau_sil = 1 in the simulation.

fit_info message Both actual and predicted relative reductions in the sum of squares are at most 0.000000

fit_info param_cov from scipy [[ 7.71439959e-03 -1.13153353e-04 -3.08345771e-03 -9.72135409e-06 8.85727485e-04 5.37435623e-03] [-1.13153353e-04 3.95494604e-04 -5.27363936e-05 -3.62103846e-06 -5.54130404e-06 -3.62487661e-04] [-3.08345771e-03 -5.27363936e-05 5.54494824e-03 1.72732595e-04 -3.66794776e-03 -3.67701879e-03] [-9.72135409e-06 -3.62103846e-06 1.72732595e-04 2.21930261e-05 -3.87891140e-04 -4.76215557e-04] [ 8.85727485e-04 -5.54130404e-06 -3.66794776e-03 -3.87891140e-04 7.77305065e-03 1.17720169e-02] [ 5.37435623e-03 -3.62487661e-04 -3.67701879e-03 -4.76215557e-04 1.17720169e-02 3.13437185e-02]]

astropy modeling uncertainties standard deviations amplitude_0| 0.088 x_0_0 | 0.02 fwhm_0 | 0.074 slope_1 | 0.005 intercept_1| 0.088 tau_sil_2 | 0.177

astropy modeling parameters Model: CompoundModel Inputs: ('x',) Outputs: ('y',) Model set size: 1 Expression: ([0] + [1]) * [2] Components: [0]: <Drude1D(amplitude=3.2482878, x_0=6.24199046, fwhm=1.5904412)>

[1]: <Linear1D(slope=0.00508949, intercept=0.85773984)>

[2]: <S07_attenuation(tau_sil=0.97827485)>

Parameters: amplitude_0 x_0_0 fwhm_0 slope_1 intercept_1 tau_sil_2

---

3.2482877966111023 6.241990455819625 1.5904412008494448 0.005089491650168936 0.8577398367958 0.9782748466502099

Example 2 of fitting a Compound model with (Drude + line) * S07_attenuation model (see code below). The tau_sil = 0 in the simulation.

fit_info message Both actual and predicted relative reductions in the sum of squares are at most 0.000000

fit_info param_cov from scipy None

astropy modeling uncertainties None

astropy modeling parameters Model: CompoundModel Inputs: ('x',) Outputs: ('y',) Model set size: 1 Expression: ([0] + [1]) * [2] Components: [0]: <Drude1D(amplitude=3.21581905, x_0=6.24556186, fwhm=1.51565112)>

[1]: <Linear1D(slope=-0.00201203, intercept=0.97010531)>

[2]: <S07_attenuation(tau_sil=0.)>

Parameters: amplitude_0 x_0_0 fwhm_0 slope_1 intercept_1 tau_sil_2

---

3.215819053682457 6.245561860226392 1.5156511177381597 -0.0020120306841257284 0.9701053054797406       0.0

Test code used (based on the code given in PR #143):

import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling import models, fitting
from astropy.modeling.physical_models import Drude1D, BlackBody
from astropy.modeling import Fittable1DModel
from astropy.modeling import Parameter
from astropy.modeling.models import custom_model

from pahfit.component_models import S07_attenuation

import astropy.units as u

def main(prof='Gauss', plot=False):
    # Generate fake data
    np.random.seed(0)
    x = np.linspace(1.0, 20., 200)
    y = 3 * np.exp(-0.5 * (x - 6.3)**2 / 0.8**2)
    y += 1.0   # and offset for line to fit
    y += np.random.normal(0., 0.2, x.shape)
    y_unc = np.full((len(y)), 0.2)

    # attenuation
    att = S07_attenuation(tau_sil=1.0)

    # y *= att(x)   # multiple by the dust attenuation

    if prof == 'Gauss':
        # Fit the data using a Gaussian profile
        g_init = models.Gaussian1D(amplitude=1., mean=6.0, stddev=1.)
        fit_g = fitting.LevMarLSQFitter(calc_uncertainties=True)
        g = fit_g(g_init, x, y)
    elif prof == 'Drude':
        # Fit the data using a Drude profile
        g_init = models.Drude1D(amplitude=3., x_0=6.3, fwhm=1.)
        fit_g = fitting.LevMarLSQFitter(calc_uncertainties=True)
        g = fit_g(g_init, x, y)

    # Define line model
    line = models.Linear1D(slope=0.0001)

    # Combine models

    # Gaussian + line
    comp1 = (g_init + line) * att
    fit_m1 = fitting.LevMarLSQFitter(calc_uncertainties=True)
    comp_fit1 = fit_m1(comp1, x, y, weights=1.0 / y_unc)

    print("fit_info message")
    print(fit_m1.fit_info["message"])

    print("fit_info param_cov from scipy")
    print(fit_m1.fit_info["param_cov"])

    print("astropy modeling uncertainties")
    print(comp_fit1.stds)

    print("astropy modeling parameters")
    print(comp_fit1)

    if plot:
        # Plot the data with the best-fit model
        plt.figure(figsize=(8, 5))
        plt.plot(x, y, 'ko')
        plt.plot(x, g_init(x), label=f'{prof} init')
        plt.plot(x, g(x), label=f'{prof}')
        plt.plot(x, comp_fit1(x), label=f'{prof}+line')
        plt.xlabel('Position')
        plt.ylabel('Flux')
        plt.legend(loc=2)

        plt.show()

if __name__ == '__main__':
    main(prof='Drude', plot=True)

[karllark]

I'm not quite sure how to proceed. Do others have more knowledge of the LevMar algorithm that could provide insight? Maybe there are parameters we could set? @jdtsmith?

One possible longer term solution would be to use the LevMar fit as the starting point for a MCMC sampler (e.g., emcee) and then compute the parameter uncertainties (and any combined parameter uncertainties) needed from the sampler chains. This should be more robust and provide full 1d posterior probability distribution functions (pPDFs) as well. Or 2D pPDFs. Or whatever nD pPDFs one is interested in.

[karllark]

Current astropy.modeling LevMarLSQFitter uses older scipy.optimize.leastsq interface. https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

If updating to new scipy.optimize.least_squares interface would help (e.g., access to more tuning parameters from what I can see), this could be done (I could put in a PR to astropy to do this). https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

Both use the MINPACK implementation of the Levenberg-Marquardt algorithm.

[jdtsmith]

LM certainly provides uncertainty estimates as a matter of course, since it uses these during optimization to make informed guesses about steps to take. BUT, I know that MPFIT for example sets parameter uncertainty to zero when it hits a boundary:

PERROR - The formal 1-sigma errors in each parameter, computed 
            from the covariance matrix. If a parameter is held 
            fixed, or if it touches a boundary, then the error is 
            reported as zero. 

You can see why they might make that choice: if you've hit a boundary, by definition you now have a limiting result. I believe @ThomasSYLai asked me that questions about non-zero parameters in the IDL PAHFIT with zero uncertainty on Monday, and I think this is the answer. If scipy makes the same choice, this could be the problem. What happens if you remove all fitting parameter limits?

[alexmaragko]

Current astropy.modeling LevMarLSQFitter uses older scipy.optimize.leastsq interface. https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

If updating to new scipy.optimize.least_squares interface would help (e.g., access to more tuning parameters from what I can see), this could be done (I could put in a PR to astropy to do this). https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

Both use the MINPACK implementation of the Levenberg-Marquardt algorithm.

We can check if for example scipy's curve_fit would return the covariance matrix in our test examples. Setting something like: pars, covar = curve_fit(func, xdata, ydata, sigma=ydata_unc, p0=p0, bounds=bounds, absolute_sigma=True) where the func would be our compound model. And then retrieve uncertainties by: perr = np.sqrt(np.diag(covar)).