13.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.

13.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.

13.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 13.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.

13.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.