diff --git a/examples/doc_builtinmodels_nistgauss.py b/examples/doc_builtinmodels_nistgauss.py
index 67a5e4e4..c206c754 100644
--- a/examples/doc_builtinmodels_nistgauss.py
+++ b/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)
diff --git a/examples/doc_builtinmodels_peakmodels.py b/examples/doc_builtinmodels_peakmodels.py
index 0708a252..7e4a935d 100644
--- a/examples/doc_builtinmodels_peakmodels.py
+++ b/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))
diff --git a/examples/doc_confidence_chi2_maps.py b/examples/doc_confidence_chi2_maps.py
index 4f69b81e..b537fda3 100644
--- a/examples/doc_confidence_chi2_maps.py
+++ b/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:
diff --git a/examples/example_fit_with_bounds.py b/examples/example_fit_with_bounds.py
index 2e176b25..96453c7a 100644
--- a/examples/example_fit_with_bounds.py
+++ b/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')
diff --git a/lmfit/confidence.py b/lmfit/confidence.py
index bac0490a..4396fc1a 100644
--- a/lmfit/confidence.py
+++ b/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
diff --git a/lmfit/minimizer.py b/lmfit/minimizer.py
index b45ae8c1..4ac74ce5 100644
--- a/lmfit/minimizer.py
+++ b/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)
diff --git a/lmfit/model.py b/lmfit/model.py
index 9c30b5ef..9d36166d 100644
--- a/lmfit/model.py
+++ b/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,6 +1663,11 @@ 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
-------
@@ -1670,9 +1675,10 @@ def fit_report(self, modelpars=None, show_correl=True,
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}'
diff --git a/lmfit/printfuncs.py b/lmfit/printfuncs.py
index d77e4ab4..dbdc98ba 100644
--- a/lmfit/printfuncs.py
+++ b/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'Correlations {extra}
')
add('')
for name1, name2, val in sort_correls:
- stat_row(name1, name2, f"{val:.4f}")
+ stat_row(name1, name2, f"{val:+.4f}")
add('
')
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.
diff --git a/tests/test_covariance_matrix.py b/tests/test_covariance_matrix.py
index fbf92354..8c7d7336 100644
--- a/tests/test_covariance_matrix.py
+++ b/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)
diff --git a/tests/test_nose.py b/tests/test_nose.py
index 7980d4a1..b8237f5f 100644
--- a/tests/test_nose.py
+++ b/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):
diff --git a/tests/test_printfuncs.py b/tests/test_printfuncs.py
index 7b8f1341..30314e7d 100644
--- a/tests/test_printfuncs.py
+++ b/tests/test_printfuncs.py
@@ -7,8 +7,8 @@
report_ci, report_fit)
from lmfit.lineshapes import gaussian
from lmfit.models import GaussianModel
-from lmfit.printfuncs import (alphanumeric_sort, fitreport_html_table,
- getfloat_attr, gformat)
+from lmfit.printfuncs import (alphanumeric_sort, correl_table,
+ fitreport_html_table, getfloat_attr, gformat)
np.random.seed(0)
@@ -202,6 +202,22 @@ def sort_length(s):
assert 'fwhm' in report_split[indx_vars+1]
+def test_correl_table(fitresult, capsys):
+ """Verify that ``correl_table`` is not empty."""
+ table_lines = correl_table(fitresult.params).split('\n')
+ nvarys = fitresult.nvarys
+
+ assert len(table_lines) == nvarys+4
+ assert len(table_lines[5]) > nvarys*10
+
+
+def test_fit_report_correl_table(fitresult, capsys):
+ """Verify that ``correl_table`` is not empty."""
+ out = fitresult.fit_report(correl_mode='table')
+ assert '[[Correlations]]' in out
+ assert '----+' in out
+
+
def test_report_fit(fitresult, capsys):
"""Verify that the fit report is printed when using report_fit."""
# report_fit with MinimizerResult/ModelResult as argument gives full