Merge pull request #831 from lmfit/calc_dely_components

add calculation of dely for components of a composite model

azure-pipelines.yml

@@ -71,8 +71,6 @@ stages:
             python.version: '3.7'
           Python38:
             python.version: '3.8'
-          Python39:
-            python.version: '3.9'
 
       steps:
         - task: UsePythonVersion@0
@@ -86,7 +84,7 @@ stages:
           displayName: 'Install dependencies'
         - script: |
             python -m pip install --upgrade build pip wheel
-            python -m pip install asteval==0.9.22 numpy==1.18 scipy==1.4.0 uncertainties==3.1.4
+            python -m pip install asteval==0.9.22 numpy==1.18 scipy==1.4.0 uncertainties==3.1.4 pytest coverage codecov flaky
           displayName: 'Install minimum required version of dependencies'
         - script: |
             python -m build

doc/model.rst

@@ -714,6 +714,17 @@ and ``bic``.
 
    numpy.ndarray of data to compare to model.
 
+.. attribute:: dely
+
+   numpy.ndarray of estimated uncertainties in the ``y`` values of the model
+   from :meth:`ModelResult.eval_uncertainty`  (see :ref:`eval_uncertainty_sec`).
+
+.. attribute:: dely_comps
+
+   a dictionary of estimated uncertainties in the ``y`` values of the model
+   components, from :meth:`ModelResult.eval_uncertainty` (see
+   :ref:`eval_uncertainty_sec`).
+
 .. attribute:: errorbars
 
    Boolean for whether error bars were estimated by fit.
@@ -822,6 +833,8 @@ and ``bic``.
    array, so that ``weights*(data - fit)`` is minimized in the
    least-squares sense.
 
+.. _eval_uncertainty_sec:
+
 Calculating uncertainties in the model function
 -----------------------------------------------
 
@@ -861,6 +874,21 @@ figure below.
     plt.show()
 
 
+.. versionadded:: 1.0.4
+
+If the model is a composite built from multiple components, the
+:meth:`ModelResult.eval_uncertainty` method will evaluate the uncertainty of
+both the full model (often the sum of multiple components) as well as the
+uncertainty in each component.  The uncertainty of the full model will be held in
+``result.dely``, and the uncertainties for each component will be held in the dictionary
+``result.dely_comps``, with keys that are the component prefixes.
+
+An example script shows how the uncertainties in components of a composite
+model can be calculated and used:
+
+.. jupyter-execute:: ../examples/doc_model_uncertainty2.py
+
+
 .. _modelresult_saveload_sec:
 
 Saving and Loading ModelResults

doc/whatsnew.rst

@@ -16,6 +16,14 @@ consult the `lmfit GitHub repository`_.
 Version 1.0.4 Release Notes (unreleased)
 ========================================
 
+New features:
+
+- add calculation of ``dely`` for model components of composite models (Issue #761; PR #826)
+- add R^2 ``rsquared`` statistic to fit outputs and reports for Model fits (Issue #803; PR #810)
+- add ``SplineModel`` (PR #804)
+- add ``Pearson4Model`` (@lellid; PR #800)
+
+
 Bug fixes/enhancements:
 
 - make sure variable ``spercent`` is always defined in ``params_html_table`` functions (reported by @MySlientWind; Issue #768, PR #770)
@@ -24,8 +32,10 @@ Bug fixes/enhancements:
 - components used to create a ``CompositeModel`` can now have different independent variables (@Julian-Hochhaus; Discussion #787; PR #788))
 - fixed function definition for ``StepModel(form='linear')``, was not consistent with the other ones. (@matpompili; PR #794)
 - fixed height factor for ``Gaussian2dModel``, was not correct. (@matpompili; PR #795)
-- added ``Pearson4`` fitting model
 - for covariances with negative diagonal elements, we set the covariance to ``None`` (PR #813)
+- fixed linear mode for ``RectangleModel`` (@arunpersaud; Issue #815; PR #816)
+- report correct initial values for parameters with bounds (Issue #820; PR #821)
+- allow recalculation of confidence intervals (@jagerber48; PR #798)
 
 Deprecations:
 

examples/doc_model_uncertainty2.py

@@ -0,0 +1,84 @@
+# <examples/doc_model_uncertainty2.py>
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lmfit.models import ExponentialModel, GaussianModel
+
+dat = np.loadtxt('NIST_Gauss2.dat')
+x = dat[:, 1]
+y = dat[:, 0]
+
+model = (GaussianModel(prefix='g1_') +
+         GaussianModel(prefix='g2_') +
+         ExponentialModel(prefix='bkg_'))
+
+params = model.make_params(bkg_amplitude=100, bkg_decay=80,
+                           g1_amplitude=3000,
+                           g1_center=100,
+                           g1_sigma=10,
+                           g2_amplitude=3000,
+                           g2_center=150,
+                           g2_sigma=10)
+
+result = model.fit(y, params, x=x)
+
+print(result.fit_report(min_correl=0.5))
+
+
+comps = result.eval_components(x=x)
+dely = result.eval_uncertainty(sigma=3)
+
+
+fig, axes = plt.subplots(2, 2, figsize=(12.8, 9.6))
+
+axes[0][0].plot(x, y, 'o', color='#99002299', markersize=3, label='data')
+axes[0][0].plot(x, result.best_fit, '-', label='best fit')
+axes[0][0].plot(x, result.init_fit, '--', label='initial fit')
+axes[0][0].set_title('data, initial fit, and best-fit')
+axes[0][0].legend()
+
+axes[0][1].plot(x, y, 'o', color='#99002299', markersize=3, label='data')
+axes[0][1].plot(x, result.best_fit, '-', label='best fit')
+axes[0][1].fill_between(x, result.best_fit-dely, result.best_fit+dely,
+                        color="#8A8A8A", label=r'3-$\sigma$ band')
+axes[0][1].set_title('data, best-fit, and uncertainty band')
+axes[0][1].legend()
+
+
+axes[1][0].plot(x, result.best_fit, '-', label=r'best fit, 3-$\sigma$ band')
+axes[1][0].fill_between(x,
+                        result.best_fit-result.dely,
+                        result.best_fit+result.dely,
+                        color="#8A8A8A")
+
+axes[1][0].plot(x, comps['bkg_'], label=r'background, 3-$\sigma$ band')
+axes[1][0].fill_between(x,
+                        comps['bkg_']-result.dely_comps['bkg_'],
+                        comps['bkg_']+result.dely_comps['bkg_'],
+                        color="#8A8A8A")
+
+axes[1][0].plot(x, comps['g1_'], label=r'Gaussian #1, 3-$\sigma$ band')
+axes[1][0].fill_between(x,
+                        comps['g1_']-result.dely_comps['g1_'],
+                        comps['g1_']+result.dely_comps['g1_'],
+                        color="#8A8A8A")
+
+axes[1][0].plot(x, comps['g2_'], label=r'Gaussian #2, 3-$\sigma$ band')
+axes[1][0].fill_between(x,
+                        comps['g2_']-result.dely_comps['g2_'],
+                        comps['g2_']+result.dely_comps['g2_'],
+                        color="#8A8A8A")
+axes[1][0].set_title('model components with uncertainty bands')
+axes[1][0].legend()
+#
+axes[1][1].plot(x, result.best_fit, '-', label='best fit')
+axes[1][1].plot(x, 10*result.dely, label=r'3-$\sigma$ total (x10)')
+axes[1][1].plot(x, 10*result.dely_comps['bkg_'], label=r'3-$\sigma$ background (x10)')
+axes[1][1].plot(x, 10*result.dely_comps['g1_'], label=r'3-$\sigma$ Gaussian #1 (x10)')
+axes[1][1].plot(x, 10*result.dely_comps['g2_'], label=r'3-$\sigma$ Gaussian #2 (x10)')
+axes[1][1].set_title('uncertainties for model components')
+axes[1][1].legend()
+
+plt.show()
+
+# <end examples/doc_model_uncertainty2.py>

lmfit/model.py

@@ -7,6 +7,7 @@
 import operator
 import warnings
 
+from asteval import valid_symbol_name
 import numpy as np
 from scipy.special import erf
 from scipy.stats import t
@@ -265,7 +266,12 @@ def __init__(self, func, independent_vars=None, param_names=None,
 
         """
         self.func = func
+        if not isinstance(prefix, str):
+            prefix = ''
+        if len(prefix) > 0 and not valid_symbol_name(prefix):
+            raise ValueError(f"'{prefix}' is not a valid Model prefix")
         self._prefix = prefix
+
         self._param_root_names = param_names  # will not include prefixes
         self.independent_vars = independent_vars
         self._func_allargs = []
@@ -1497,6 +1503,10 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
            < 1, it is interpreted as the probability itself. That is,
            ``sigma=1`` and ``sigma=0.6827`` will give the same results,
            within precision errors.
+        3. Also sets attributes of `dely` for the uncertainty of the model
+           (which will be the same as the array returned by this method) and
+           `dely_comps`, a dictionary of `dely` for each component.
+
 
         Examples
         --------
@@ -1518,36 +1528,52 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
         # ensure fjac and df2 are correct size if independent var updated by kwargs
         ndata = self.model.eval(params, **userkws).size
         covar = self.covar
-        fjac = np.zeros((nvarys, ndata))
-        df2 = np.zeros(ndata)
         if any(p.stderr is None for p in params.values()):
-            return df2
+            return np.zeros(ndata)
+
+        fjac = {'0': np.zeros((nvarys, ndata))}  # '0' signify 'Full', an invalid prefix
+        df2 = {'0': np.zeros(ndata)}
+
+        for comp in self.components:
+            label = comp.prefix if len(comp.prefix) > 1 else comp._name
+            fjac[label] = np.zeros((nvarys, ndata))
+            df2[label] = np.zeros(ndata)
 
         # find derivative by hand!
         pars = params.copy()
         for i in range(nvarys):
             pname = self.var_names[i]
             val0 = pars[pname].value
             dval = pars[pname].stderr/3.0
-
             pars[pname].value = val0 + dval
-            res1 = self.model.eval(pars, **userkws)
+            res1 = {'0': self.model.eval(pars, **userkws)}
+            res1.update(self.model.eval_components(params=pars, **userkws))
 
             pars[pname].value = val0 - dval
-            res2 = self.model.eval(pars, **userkws)
+            res2 = {'0': self.model.eval(pars, **userkws)}
+            res2.update(self.model.eval_components(params=pars, **userkws))
 
             pars[pname].value = val0
-            fjac[i] = (res1 - res2) / (2*dval)
+            for key in fjac:
+                fjac[key][i] = (res1[key] - res2[key]) / (2*dval)
 
         for i in range(nvarys):
             for j in range(nvarys):
-                df2 += fjac[i]*fjac[j]*covar[i, j]
+                for key in fjac:
+                    df2[key] += fjac[key][i] * fjac[key][j] * covar[i, j]
 
         if sigma < 1.0:
             prob = sigma
         else:
             prob = erf(sigma/np.sqrt(2))
-        return np.sqrt(df2) * t.ppf((prob+1)/2.0, self.ndata-nvarys)
+
+        scale = t.ppf((prob+1)/2.0, self.ndata-nvarys)
+        self.dely = scale * np.sqrt(df2.pop('0'))
+
+        self.dely_comps = {}
+        for key in df2:
+            self.dely_comps[key] = scale * np.sqrt(df2[key])
+        return self.dely
 
     def conf_interval(self, **kwargs):
         """Calculate the confidence intervals for the variable parameters.
@@ -1656,7 +1682,7 @@ def dumps(self, **kws):
                      'ci_out', 'col_deriv', 'covar', 'errorbars', 'flatchain',
                      'ier', 'init_values', 'lmdif_message', 'message',
                      'method', 'nan_policy', 'ndata', 'nfev', 'nfree',
-                     'nvarys', 'redchi', 'residual', 'rsquared', 'scale_covar',
+                     'nvarys', 'redchi', 'rsquared', 'scale_covar',
                      'calc_covar', 'success', 'userargs', 'userkws', 'values',
                      'var_names', 'weights', 'user_options'):
 

tests/test_model_saveload.py

@@ -144,7 +144,7 @@ def test_save_load_modelresult(dill):
     text = ''
     with open(SAVE_MODELRESULT) as fh:
         text = fh.read()
-    assert 18000 < len(text) < 21000  # depending on whether dill is present
+    assert 12000 < len(text) < 60000  # depending on whether dill is present
 
     # load the saved ModelResult from file and compare results
     result_saved = load_modelresult(SAVE_MODELRESULT)

tests/test_model_uncertainties.py

@@ -1,10 +1,11 @@
 """Tests of ModelResult.eval_uncertainty()"""
+import os
 
 import numpy as np
 from numpy.testing import assert_allclose
 
 from lmfit.lineshapes import gaussian
-from lmfit.models import GaussianModel, LinearModel
+from lmfit.models import ExponentialModel, GaussianModel, LinearModel
 
 
 def get_linearmodel(slope=0.8, intercept=0.5, noise=1.5):
@@ -93,3 +94,36 @@ def test_gauss_noiselevel():
     dely_hinoise = ret2.eval_uncertainty(sigma=1)
 
     assert_allclose(dely_hinoise.mean(), 10*dely_lonoise.mean(), rtol=1.e-2)
+
+
+def test_component_uncertainties():
+    "test dely_comps"
+    y, x = np.loadtxt(os.path.join(os.path.dirname(__file__), '..',
+                                   'examples', 'NIST_Gauss2.dat')).T
+    model = (GaussianModel(prefix='g1_') +
+             GaussianModel(prefix='g2_') +
+             ExponentialModel(prefix='bkg_'))
+
+    params = model.make_params(bkg_amplitude=100, bkg_decay=80,
+                               g1_amplitude=3000,
+                               g1_center=100,
+                               g1_sigma=10,
+                               g2_amplitude=3000,
+                               g2_center=150,
+                               g2_sigma=10)
+
+    result = model.fit(y, params, x=x)
+    comps = result.eval_components(x=x)
+    dely = result.eval_uncertainty(sigma=3)
+
+    assert 'g1_' in comps
+    assert 'g2_' in comps
+    assert 'bkg_' in comps
+    assert dely.mean() > 0.8
+    assert dely.mean() < 2.0
+    assert result.dely_comps['g1_'].mean() > 0.5
+    assert result.dely_comps['g1_'].mean() < 1.5
+    assert result.dely_comps['g2_'].mean() > 0.5
+    assert result.dely_comps['g2_'].mean() < 1.5
+    assert result.dely_comps['bkg_'].mean() > 0.5
+    assert result.dely_comps['bkg_'].mean() < 1.5