Merge pull request #856 from lmfit/more_chi2_reporting_tweaks

More chi2 reporting tweaks

examples/doc_builtinmodels_nistgauss.py

@@ -26,7 +26,7 @@
 init = mod.eval(pars, x=x)
 out = mod.fit(y, pars, x=x)
 
-print(out.fit_report(min_correl=0.5))
+print(out.fit_report(correl_mode='table'))
 
 fig, axes = plt.subplots(1, 2, figsize=(12.8, 4.8))
 axes[0].plot(x, y)

examples/doc_builtinmodels_peakmodels.py

@@ -14,7 +14,7 @@
 pars = mod.guess(y, x=x)
 out = mod.fit(y, pars, x=x)
 
-print(out.fit_report(min_correl=0.25))
+print(out.fit_report(correl_mode='table'))
 
 plt.plot(x, y)
 plt.plot(x, out.best_fit, '-', label='Gaussian Model')
@@ -27,7 +27,7 @@
 pars = mod.guess(y, x=x)
 out = mod.fit(y, pars, x=x)
 
-print(out.fit_report(min_correl=0.25))
+print(out.fit_report(correl_mode='table'))
 
 plt.figure()
 plt.plot(x, y, '-')
@@ -41,7 +41,7 @@
 pars = mod.guess(y, x=x)
 out = mod.fit(y, pars, x=x)
 
-print(out.fit_report(min_correl=0.25))
+print(out.fit_report(correl_mode='table'))
 
 fig, axes = plt.subplots(1, 2, figsize=(12.8, 4.8))
 

examples/doc_confidence_chi2_maps.py

@@ -27,7 +27,7 @@
                          slope=0, intecept=2)
 
 out = mod.fit(y, params, x=x)
-print(out.fit_report())
+print(out.fit_report(correl_mode='table'))
 
 #########################
 # run conf_intervale, print report
@@ -51,24 +51,29 @@
 
 aix, aiy = 0, 0
 nsamples = 50
+explicitly_calculate_sigma = True
 
 for pairs in (('sigma', 'amplitude'), ('intercept', 'amplitude'),
               ('slope', 'intercept'), ('slope', 'center'), ('sigma', 'center')):
 
     xpar, ypar = pairs
-    print("Generating chi-square map for ", pairs)
-    c_x, c_y, dchi2_mat = conf_interval2d(out, out, xpar, ypar,
-                                          nsamples, nsamples, nsigma=3.5,
-                                          chi2_out=True)
-
-    # sigma matrix: sigma increases chi_square
-    # from  chi_square_best
-    # to    chi_square + sigma**2 * reduced_chi_square
-    # so:   sigma = sqrt(dchi2 / reduced_chi_square)
-    sigma_mat = np.sqrt(abs(dchi2_mat)/out.redchi)
-
-    # you could calculate the matrix of probabilities from sigma as:
-    # prob_mat  = np.erf(sigma_mat/np.sqrt(2))
+    if explicitly_calculate_sigma:
+        print("Generating chi-square map for ", pairs)
+        c_x, c_y, chi2_mat = conf_interval2d(out, out, xpar, ypar,
+                                             nsamples, nsamples, nsigma=3.5,
+                                             chi2_out=True)
+        # explicitly calculate sigma matrix: sigma increases chi_square
+        # from  chi_square_best
+        # to    chi_square + sigma**2 * reduced_chi_square
+        # so:   sigma = sqrt((chi2-chi2_best)/ reduced_chi_square)
+        chi2_min = chi2_mat.min()
+        sigma_mat = np.sqrt((chi2_mat-chi2_min)/out.redchi)
+    else:
+        print("Generating sigma map for ", pairs)
+        # or, you could just calculate the matrix of probabilities as:
+        # print("Generating chi-square map for ", pairs)
+        c_x, c_y, sigma_mat = conf_interval2d(out, out, xpar, ypar,
+                                              nsamples, nsamples, nsigma=3.5)
 
     aix += 1
     if aix == 2:

examples/example_fit_with_bounds.py

@@ -51,7 +51,7 @@ def residual(pars, x, data=None):
 fit = residual(out.params, x)
 
 ###############################################################################
-report_fit(out, modelpars=p_true)
+report_fit(out, modelpars=p_true, correl_mode='table')
 
 ###############################################################################
 plt.plot(x, data, 'o', label='data')

lmfit/confidence.py

@@ -370,20 +370,18 @@ def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10,
     y_name : str
         The name of the parameter which will be the y direction.
     nx : int, optional
-        Number of points in the x direction. [default = 10]
+        Number of points in the x direction (default is 10).
     ny : int, optional
-        Number of points in the y direction. [default = 10]
+        Number of points in the y direction (default is 10).
     limits : tuple, optional
         Should have the form ``((x_upper, x_lower), (y_upper, y_lower))``.
         If not given, the default is nsigma*stderr in each direction.
-    prob_func : None or callable, optional
-        Function to calculate the probability from the optimized chi-square.
-        Default is None and uses the built-in function `f_compare`
-        (i.e., F-test).
+    prob_func : None or callable, deprecated
+        Starting with version 1.2, this argument is unused and has no effect.
     nsigma : float or int, optional
-        multiplier of stderr for limits.  [default = 5.0]
+        Multiplier of stderr for limits (default is 5).
     chi2_out: bool
-        whether to return chi-square at each coordinate instead of probability.
+        Whether to return chi-square at each coordinate instead of probability.
 
     Returns
     -------
@@ -393,7 +391,7 @@ def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10,
         Y-coordinates (same shape as `ny`).
     grid : numpy.ndarray
         2-D array (with shape ``(nx, ny)``) containing the calculated
-        probabilities or chi-square
+        probabilities or chi-square.
 
     See Also
     --------
@@ -407,17 +405,16 @@ def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10,
     >>> plt.contour(x,y,gr)
 
     """
+    if prob_func is not None:
+        msg = "'prob_func' has no effect and will be removed in version 1.4."
+        raise DeprecationWarning(msg)
+
     params = result.params
 
-    best_chi = result.chisqr
+    best_chisqr = result.chisqr
+    redchi = result.redchi
     org = copy_vals(result.params)
 
-    def chi2_compare(best, current):
-        return current.chisqr - best.chisqr
-
-    if prob_func is None:
-        prob_func = chi2_compare if chi2_out else f_compare
-
     x = params[x_name]
     y = params[y_name]
 
@@ -432,29 +429,30 @@ def chi2_compare(best, current):
     y_points = np.linspace(y_lower, y_upper, ny)
     grid = np.dstack(np.meshgrid(x_points, y_points))
 
-    x.vary = False
-    y.vary = False
+    x.vary, y.vary = False, False
 
-    def calc_prob(vals, restore=False):
-        """Calculate the probability."""
-        if restore:
-            restore_vals(org, result.params)
-        x.value = vals[0]
-        y.value = vals[1]
-        save_x = result.params[x.name]
-        save_y = result.params[y.name]
-        result.params[x.name] = x
-        result.params[y.name] = y
+    def calc_chisqr(vals, restore=False):
+        """Calculate chi-square for a set of parameter values."""
+        save_x = x.value
+        save_y = y.value
+        result.params[x.name].value = vals[0]
+        result.params[y.name].value = vals[1]
         minimizer.prepare_fit(params=result.params)
         out = minimizer.leastsq()
-        prob = prob_func(result, out)
-        result.params[x.name] = save_x
-        result.params[y.name] = save_y
-        return prob
+        result.params[x.name].value = save_x
+        result.params[y.name].value = save_y
+        return out.chisqr
 
-    out = x_points, y_points, np.apply_along_axis(calc_prob, -1, grid)
+    # grid of chi-square
+    out_mat = np.apply_along_axis(calc_chisqr, -1, grid)
+
+    # compute grid of sigma values from chi-square
+    if not chi2_out:
+        chisqr0 = out_mat.min()
+        chisqr0 = min(best_chisqr, chisqr0)
+        out_mat = np.sqrt((out_mat-chisqr0)/redchi)
 
     x.vary, y.vary = True, True
     restore_vals(org, result.params)
-    result.chisqr = best_chi
-    return out
+    result.chisqr = best_chisqr
+    return x_points, y_points, out_mat

lmfit/minimizer.py

@@ -104,6 +104,10 @@ def asteval_with_uncertainties(*vals, **kwargs):
         return 0
     for val, name in zip(vals, _names):
         _asteval.symtable[name] = val
+
+    # re-evaluate all constraint parameters to
+    # force the propagation of uncertainties
+    [p._getval() for p in _pars.values()]
     return _asteval.eval(_obj._expr_ast)
 
 

lmfit/model.py

@@ -1642,7 +1642,7 @@ def ci_report(self, with_offset=True, ndigits=5, **kwargs):
                          with_offset=with_offset, ndigits=ndigits)
 
     def fit_report(self, modelpars=None, show_correl=True,
-                   min_correl=0.1, sort_pars=False):
+                   min_correl=0.1, sort_pars=False, correl_mode='list'):
         """Return a printable fit report.
 
         The report contains fit statistics and best-fit values with
@@ -1663,16 +1663,22 @@ def fit_report(self, modelpars=None, show_correl=True,
             listed in the order as they were added to the Parameters
             dictionary. If callable, then this (one argument) function is
             used to extract a comparison key from each list element.
+        correl_mode : {'list', table'} str, optional
+            Mode for how to show correlations. Can be either 'list' (default)
+            to show a sorted (if ``sort_pars`` is True) list of correlation
+            values, or 'table' to show a complete, formatted table of
+            correlations.
 
         Returns
         -------
         str
             Multi-line text of fit report.
 
         """
-        report = fit_report(self, modelpars=modelpars,
-                            show_correl=show_correl,
-                            min_correl=min_correl, sort_pars=sort_pars)
+        report = fit_report(self, modelpars=modelpars, show_correl=show_correl,
+                            min_correl=min_correl, sort_pars=sort_pars,
+                            correl_mode=correl_mode)
+
         modname = self.model._reprstring(long=True)
         return f'[[Model]]\n    {modname}\n{report}'
 

lmfit/printfuncs.py

@@ -82,7 +82,7 @@ def gformat(val, length=11):
 
 
 def fit_report(inpars, modelpars=None, show_correl=True, min_correl=0.1,
-               sort_pars=False):
+               sort_pars=False, correl_mode='list'):
     """Generate a report of the fitting results.
 
     The report contains the best-fit values for the parameters and their
@@ -104,6 +104,11 @@ def fit_report(inpars, modelpars=None, show_correl=True, min_correl=0.1,
         they were added to the Parameters dictionary. If callable, then
         this (one argument) function is used to extract a comparison key
         from each list element.
+    correl_mode : {'list', table'} str, optional
+        Mode for how to show correlations. Can be either 'list' (default)
+        to show a sorted (if ``sort_pars`` is True) list of correlation
+        values, or 'table' to show a complete, formatted table of
+        correlations.
 
     Returns
     -------
@@ -191,7 +196,11 @@ def fit_report(inpars, modelpars=None, show_correl=True, min_correl=0.1,
         else:
             add(f"    {nout} {par.value: .7g} (fixed)")
 
-    if show_correl:
+    if show_correl and correl_mode.startswith('tab'):
+        add('[[Correlations]] ')
+        for line in correl_table(params).split('\n'):
+            buff.append('  %s' % line)
+    elif show_correl:
         correls = {}
         for i, name in enumerate(parnames):
             par = params[name]
@@ -211,7 +220,7 @@ def fit_report(inpars, modelpars=None, show_correl=True, min_correl=0.1,
             maxlen = max(len(k) for k in list(correls.keys()))
         for name, val in sort_correl:
             lspace = max(0, maxlen - len(name))
-            add(f"    C({name}){(' '*30)[:lspace]} = {val:.3f}")
+            add(f"    C({name}){(' '*30)[:lspace]} = {val:+.4f}")
     return '\n'.join(buff)
 
 
@@ -274,11 +283,47 @@ def stat_row(label, val, val2=''):
             add(f'<h2>Correlations {extra}</h2>')
             add('<table>')
             for name1, name2, val in sort_correls:
-                stat_row(name1, name2, f"{val:.4f}")
+                stat_row(name1, name2, f"{val:+.4f}")
             add('</table>')
     return ''.join(html)
 
 
+def correl_table(params):
+    """Return a printable correlation table for a Parameters object."""
+    varnames = [vname for vname in params if params[vname].vary]
+    nwid = max(8, max([len(vname) for vname in varnames])) + 1
+
+    def sfmt(a):
+        return f" {a:{nwid}s}"
+
+    def ffmt(a):
+        return sfmt(f"{a:+.4f}")
+
+    title = ['', sfmt('Variable')]
+    title.extend([sfmt(vname) for vname in varnames])
+
+    title = '|'.join(title) + '|'
+    bar = [''] + ['-'*(nwid+1) for i in range(len(varnames)+1)] + ['']
+    bar = '+'.join(bar)
+
+    buff = [bar, title, bar]
+
+    for vname, par in params.items():
+        if not par.vary:
+            continue
+        line = ['', sfmt(vname)]
+        for vother in varnames:
+            if vother == vname:
+                line.append(ffmt(1))
+            elif vother in par.correl:
+                line.append(ffmt(par.correl[vother]))
+            else:
+                line.append('unknown')
+        buff.append('|'.join(line) + '|')
+    buff.append(bar)
+    return '\n'.join(buff)
+
+
 def params_html_table(params):
     """Return an HTML representation of Parameters.
 

tests/test_covariance_matrix.py

@@ -98,8 +98,8 @@ def test_bounds_expression():
     assert_almost_equal(result.params['amplitude'].stderr, 0.13861506,
                         decimal=4)
     assert_almost_equal(result.params['gamma'].stderr, 0.00368468, decimal=4)
-    assert_almost_equal(result.params['fwhm'].stderr, 0.00806917, decimal=4)
-    assert_almost_equal(result.params['height'].stderr, 0.03009459, decimal=4)
+    assert_almost_equal(result.params['fwhm'].stderr, 0.01326862028, decimal=4)
+    assert_almost_equal(result.params['height'].stderr, 0.0395990, decimal=4)
 
     assert_almost_equal(result.params['sigma'].correl['center'],
                         -4.6623973788006615e-05, decimal=4)

tests/test_nose.py

@@ -12,6 +12,7 @@
 from lmfit.lineshapes import gaussian
 from lmfit.minimizer import (HAS_EMCEE, SCALAR_METHODS, MinimizerResult,
                              _nan_policy)
+from lmfit.models import SineModel
 
 try:
     import numdifftools  # noqa: F401
@@ -317,6 +318,32 @@ def test_ufloat():
     assert_allclose(y.std_dev, 0.0, rtol=1.e-7)
 
 
+def test_stderr_propagation():
+    """Test propagation of uncertainties to constraint expressions."""
+    model = SineModel()
+    params = model.make_params(amplitude=1, frequency=9.0, shift=0)
+    params.add("period", expr="1/frequency")
+    params.add("period_2a", expr="2*period")
+    params.add("period_2b", expr="2/frequency")
+    params.add("thing_1", expr="shift + frequency")
+    params.add("thing_2", expr="shift + 1/period")
+
+    np.random.seed(3)
+    xs = np.linspace(0, 1, 51)
+    ys = np.sin(7.45 * xs) + 0.01 * np.random.normal(size=xs.shape)
+
+    result = model.fit(ys, x=xs, params=params)
+
+    opars = result.params
+    assert_allclose(opars['period'].stderr, 1.1587e-04, rtol=1.e-3)
+    assert_allclose(opars['period_2b'].stderr, 2.3175e-04, rtol=1.e-3)
+    assert_allclose(opars['thing_1'].stderr, 0.0037291, rtol=1.e-3)
+
+    assert_allclose(opars['period_2a'].stderr, 2*opars['period'].stderr, rtol=1.e-5)
+    assert_allclose(opars['period_2b'].stderr, opars['period_2a'].stderr, rtol=1.e-5)
+    assert_allclose(opars['thing_1'].stderr, opars['thing_2'].stderr, rtol=1.e-5)
+
+
 class CommonMinimizerTest(unittest.TestCase):
 
     def setUp(self):