| Larch: Data Analysis Tools for X-ray Spectroscopy

# 12.8. Advanced Confidence Intervals and Chi-square maps

Having good estimates for uncertainties in fitted parameters is important for any scientific analysis. One of principle attractions to using the Levenberg-Marquardt algorithm that is the default fitting mechanism with the `minimize()` function is that it will automatically calculate estimates of parameter uncertainties and correlations. This is very convenient and generally reliable, but it should be made clear that the basic assumptions made when these uncertainties are estimated are not always perfect. Unfortunately, it is sometimes difficult to tell when this is the case.

It is therefore fairly common to see analyses that include explicit exploration of Parameter values away from their best-fit solution, in order to determine the degree of confidence in the best-fit values. Larch provides two main functions to help explore such cases. To be sure, they are much slower than the automatic estimation of the uncertainties. For many (perhaps most) cases, they do not provide much better insight than the automatic method.

confidence_intervals(minimizer, sigmas=(1, 2, 3), prob_func=None, maxiter=200)

  Calculate confidence intervals for the parameters from a given fit.

  Parameters

    - minimizer – the minimizer object returned by `minimize()`.

    - sigmas – list of sigma-levels to find parameter values for.

    - prob_func – `None` or callable function to calculate the probability from the optimized chi-square. By default, `f_compare()`, the standard F-test, is used.

  Returns

    a dictionary of parameter names, with each value containing a list of (sigma, value) pairs.

This function will adjust the value for each parameter, re-optimizing the other parameters until it finds the parameter values that increase sigma by the levels indicated.

confidence_report(conf_values)

  Convert the output of `confidence_intervals()` into a printable report.

  Parameters

    conf_values – confidence values returned by `confidence_intervals()`.

  Returns

    a string containing the report, which can be printed or stored.

f_test(ndata, nparams, chisquare, chisquare0, nfix=1)

  Returns the standard F-test value for the probability that one fit is better than another.

## 12.8.1. Confidence interval, Example 1

Let’s begin with a shortened version of the first example from the previous section.

    ## examples/fitting/doc_example_conf1.lar
    # create mock data
    mdat = group()
    mdat.x = linspace(-10, 10, 201)
    mdat.y = 1.0 + gaussian(mdat.x, 12.0, 1.5, 2.0) + \
             random.normal(size=len(mdat.x), scale=0.050)

    params = param_group(off = guess(0),
                         amp = guess(5, min=0),
                         cen = guess(2),
                         wid = guess(1, min=0))

    # define objective function for fit residual
    def myresid(p, data):
        return data.y - (p.off + gaussian(data.x, p.amp, p.cen, p.wid))
    #enddef

    # perform fit

    result = minimize(myresid, params, args=(mdat,))

    # print report of parameters, uncertainties
    print( fit_report(result))

    out = confidence_intervals(result, sigmas=(1, 2, 3))

    print( confidence_report(out))
    ## end of examples/fitting/doc_example_conf1.lar

The printed output from `fit_report(params)` will include this:

    [[Variables]]
       amp            =  12.129867 +/- 0.124722 (init=  5.000000)
       cen            =  1.476822 +/- 0.017266 (init=  2.000000)
       off            =  0.998814 +/- 0.007131 (init=  0.000000)
       wid            =  2.022301 +/- 0.016938 (init=  1.000000)

    [[Correlations]]     (unreported correlations are <  0.100)
       amp, off             = -0.861
       amp, wid             =  0.808
       off, wid             = -0.692
    =======================================================

while the output from the much more explicit search done in `confidence_intervals()` and reported by `confidence_report()` will be:

    # Confidence Interval Report
    # Sigmas:          -3         -2         -1          0          1          2          3
    # Percentiles:  -99.730    -95.450    -68.269      0.000     68.269     95.450     99.730
    #==========================================================================================
     amp             11.755      11.88     12.005      12.13     12.256     12.385     12.517
        -best      -0.37459   -0.24967   -0.12455          0    0.12604    0.25555    0.38759
     wid             1.9707     1.9885     2.0053     2.0223     2.0395     2.0569     2.0749
        -best     -0.051552  -0.033784  -0.016975          0   0.017156   0.034567   0.052626
     off             0.9766    0.98468     0.9912    0.99881     1.0064     1.0128     1.0209
        -best     -0.022212  -0.014138  -0.0076156          0  0.0075667   0.013987   0.022049
     cen             1.4248     1.4423     1.4596     1.4768     1.4941     1.5113     1.5289
        -best     -0.052057   -0.03452  -0.017248          0   0.017249   0.034524   0.052068

The automatic error estimates given from `minimize()` are meant to be 1-\(\sigma\) uncertainties. Comparing the two methods we find:

> Parameter | Best Value | Automatic 1-\(\sigma\) value | Explicit 1-\(\sigma\) value
> --- | --- | --- | ---
> amp | 12.1299 | +/- 0.1247 | +0.1267, -0.1246
> cen | 1.4768 | +/- 0.0173 | +0.0172, -0.0172
> off | 0.9988 | +/- 0.0071 | +0.0076, -0.0076
> cen | 2.0223 | +/- 0.0169 | +0.0172, -0.0170

which seems to justify the use of the automated method. The uncertainties found from the more thorough exploration shows symmetric uncertainties, even out to the 3-\(\sigma\) level, and of the 4 1-\(\sigma\) uncertainties, 3 are within 2%, and the worst agreement, for the smallest uncertainty is within 7%. It also shows that the scaling of uncertainties is fairly linear with \(\sigma\): the 3-\(\sigma\) values are approximately 3 times the 1-\(\sigma\) values.

## 12.8.2. Confidence interval, Example 2

Of course, there are more challenging cases than the one above. A double exponential function is one such example, so we start with a fit to mock data

    ## examples/fitting/doc_example_conf2.lar
    #  Shows usage of fitting with non-normal correlation of variables,
    #  and confidence intervals.

    random.seed(1)
    dat = group(x = linspace(0, 10, 101))
    dat.y = 3*exp(-dat.x/2) - 5*exp(-dat.x/9) + random.normal(size=len(dat.x), scale=0.05)

    fitparams = param_group(a1 = guess(3.5),
                            a2 = guess(-9.5),
                            t1 = guess(3),
                            t2 = guess(15))

    def fit_exp(pars, dat):
        model = pars.a1 * exp(-dat.x / pars.t1) + pars.a2 * exp(-dat.x / pars.t2)
        return model - dat.y
    enddef

    minout = minimize(fit_exp, fitparams, args=(dat,))
    print( fit_report(minout, min_correl=0))

    final = fit_exp(fitparams, dat) + dat.y
    newplot (dat.x, dat.y, label='data', marker='o', linewidth=0)
       plot (dat.x, final, label='fit')

    ## end of examples/fitting/doc_example_conf2.lar

The resulting statistics report with the automated uncertainties is:

    ===================== FIT RESULTS =====================
    [[Statistics]]    Fit succeeded,  method = 'leastsq'.
       Message from fit    = Fit succeeded.
       npts, nvarys, nfree = 101, 4, 97
       nfev (func calls)   = 36
       chi_square          = 0.191322
       reduced chi_square  = 0.001972

    [[Variables]]
       a1             =  2.828857 +/- 0.149776 (init=  3.500000)
       a2             = -4.819553 +/- 0.159495 (init= -9.500000)
       t1             =  1.878519 +/- 0.100212 (init=  3.000000)
       t2             =  9.270866 +/- 0.309035 (init=  15.000000)

    [[Correlations]]     (unreported correlations are <  0.100)
       a2, t2               =  0.991
       a1, a2               = -0.991
       a1, t2               = -0.988
       a2, t1               = -0.968
       t1, t2               = -0.937
       a1, t1               =  0.935
    =======================================================

You can see that the correlations between all 6 pairs of variables is above 90%. The resulting plot of the best-fit looks fairly reasonable:

Figure 12.8.2.1 Fit to double exponential function.

But if we we ask for the more thorough investigation of the confidence intervals in these parameters with:

    conf_int = confidence_intervals(minout)
    print confidence_report(conf_int)

the resulting report is:

    # Confidence Interval Report
    # Sigmas:          -3         -2         -1          0          1          2          3
    # Percentiles:  -99.730    -95.450    -68.269      0.000     68.269     95.450     99.730
    #==========================================================================================
     a1              2.4622     2.5704      2.691     2.8289     2.9936     3.1985      3.467
        -best      -0.36665   -0.25842   -0.13785          0    0.16472    0.36968    0.63812
     a2             -5.4926    -5.2111    -4.9945    -4.8196    -4.6726     -4.544     -4.429
        -best      -0.67309   -0.39152   -0.17499          0    0.14698    0.27553    0.39058
     t2              8.2604     8.6221     8.9556     9.2709     9.5761     9.8776     10.182
        -best       -1.0105   -0.64882   -0.31531          0    0.30526    0.60676    0.91132
     t1              1.6107     1.6942     1.7827     1.8785     1.9841     2.1049     2.2456
        -best      -0.26785   -0.18436   -0.09583          0    0.10553    0.22637    0.36709

Now can see more asymmetric uncertainty values, specifically that the -n-\(\sigma\) and +n-\(\sigma\) are different, and don’t seem to be linear in n. Comparing the 1-\(\sigma\) levels between the automated and explicit methods as we did above, we now have

> Parameter | Best Value | Automatic 1-\(\sigma\) value | Explicit 1-\(\sigma\) value
> --- | --- | --- | ---
> a1 | 2.8289 | +/- 0.1498 | +0.1647, -0.1379
> a2 | -4.8196 | +/- 0.1595 | +0.1470, -0.1750
> t1 | 9.2709 | +/- 0.1002 | +0.1055, -0.0958
> t2 | 1.8785 | +/- 0.3090 | +0.3053, -0.3153

In fairness, the automated values don’t look too bad, given that they cannot reflect asymmetric uncertainties. But, like the reported correlations, the full report above hints at a less than ideal case.

## 12.8.3. Chi-square Maps

We can further explore the correlation between pairs of variables by making and visualizing a map of the chi-square (\(\chi^2\)) statistic. This can be helpful to determine if the automatically calculated uncertainties and correlation are reasonable, and to look for pathological cases. The `chi2_map()` function will calculate a map of \(\chi^2\) for a pair of variable parameters by brute force.

chi2_map(minimizer, xname, yname, nx=11, ny=11, sigmas=5, xrange=None, yrange=None)

  Return a chi-square map for two parameters in a fit

  Parameters

    - minimizer – the minimizer object returned by `minimize()`.

    - xname – name of variable for x-axis

    - yname – name of variable for y-axis

    - nx – number of steps in x [11]

    - ny – number of steps in y [11]

    - sigmas – extent of x, y values to calculate, in \(\sigma\)

    - xrange – explicit x range of calculations [x.best +/- `sigmas` * x.stderr]

    - yrange – explicit y range of calculations [y.best +/- `sigmas` * y.stderr]

  Returns

    a tuple of (x_vals, y_vals, chi2_map)

  The `xrange` and `yrange` arguments can be used to fully dictate the x and y values to use. By default, the x and y values are automatically determined from `nx`, `ny`, and `sigmas`, with the `sigmas` argument sets how far from the best value to extend the ranges.

As an example usage, we return to the first example of the “well-behaved” fit, and run `chi2_map()` on a pair of variables with low correlation and a pair with high correlation:

    x1, y1, chi2_ampcen = chi2_map(mout, 'amp', 'cen', nx=21, ny=21)
    x2, y2, chi2_ampwid = chi2_map(mout, 'amp', 'wid', nx=21, ny=21)
    contour(chi2_ampcen, x=x1, y=y1, title='Correlation(amp, cen)')
    contour(chi2_ampwid, x=x2, y=y2, win=2, title='Correlation(amp, wid)')

with the resulting Chi-square maps looking like this:

The circular map for the uncorrelated parameters amp and cen and the elliptical map for the highly correlated parameters amp and wid are exactly as would be expected, and what the automated estimate of uncertainties and correlations assumes.

But now, if we turn to the more pathological case of the double exponential, we calculate the chi-square maps as:

    x1, y1, c1 = chi2_map(minout, 'a1', 'a2', nx=25, ny=25)
    x2, y2, c2 = chi2_map(minout, 'a1', 't1', nx=25, ny=25)
    contour(c1, x=x1, y=y1,        title='Correlation(a1, a2)')
    contour(c2, x=x2, y=y2, win=2, title='Correlation(a1, t1)')

with the resulting contour plots:

Here, the values of chi-square quickly grow very large away from the ideal fit. More importantly, the maps are not remotely elliptical.

Finally, it should be emphasized that while all the tests in this section are informative, they are also fairly slow, re-running the fits many times. It is probably safe to rely on the automatic calculations of uncertainties and correlations, and use these methods on occasions of extremely high correlation, or when nearing a final analysis.