Add Pearson4 fit model (original Wikipedia version)

doc/builtin_models.rst

@@ -97,6 +97,11 @@ of 0 on the value of ``sigma``.
 
 .. autoclass:: MoffatModel
 
+:class:`Pearson4Model`
+~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: Pearson4Model
+
 :class:`Pearson7Model`
 ~~~~~~~~~~~~~~~~~~~~~~
 

doc/whatsnew.rst

@@ -24,6 +24,7 @@ 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
 
 Deprecations:
 

lmfit/lineshapes.py

@@ -1,9 +1,11 @@
 """Basic model lineshapes and distribution functions."""
 
-from numpy import (arctan, copysign, cos, exp, isclose, isnan, log, pi, real,
-                   sin, sqrt, where)
+from numpy import (arctan, copysign, cos, exp, isclose, isnan, log, log1p, pi,
+                   real, sin, sqrt, where)
+from scipy.special import betaln as betalnfcn
 from scipy.special import erf, erfc
 from scipy.special import gamma as gamfcn
+from scipy.special import loggamma as loggammafcn
 from scipy.special import wofz
 
 log2 = log(2)
@@ -14,7 +16,7 @@
 tiny = 1.0e-15
 
 functions = ('gaussian', 'gaussian2d', 'lorentzian', 'voigt', 'pvoigt',
-             'moffat', 'pearson7', 'breit_wigner', 'damped_oscillator',
+             'moffat', 'pearson4', 'pearson7', 'breit_wigner', 'damped_oscillator',
              'dho', 'logistic', 'lognormal', 'students_t', 'expgaussian',
              'doniach', 'skewed_gaussian', 'skewed_voigt',
              'thermal_distribution', 'step', 'rectangle', 'exponential',
@@ -145,6 +147,26 @@ def moffat(x, amplitude=1, center=0., sigma=1, beta=1.):
     return amplitude / (((x - center)/max(tiny, sigma))**2 + 1)**beta
 
 
+def pearson4(x, amplitude=1.0, center=0.0, sigma=1.0, expon=1.0, skew=0.0):
+    """Return a Pearson4 lineshape.
+
+    Using the Wikipedia definition:
+
+    pearson4(x, amplitude, center, sigma, expon, skew) =
+        amplitude*|gamma(expon + I skew/2)/gamma(m)|**2/(w*beta(expon-0.5, 0.5)) * (1+arg**2)**(-expon) * exp(-skew * arctan(arg))
+
+    where ``arg = (x-center)/sigma``, `gamma` is the gamma function and `beta` is the beta function.
+
+    For more information, see: https://en.wikipedia.org/wiki/Pearson_distribution#The_Pearson_type_IV_distribution
+
+    """
+    expon = max(tiny, expon)
+    sigma = max(tiny, sigma)
+    arg = (x - center) / sigma
+    logprefactor = 2 * (real(loggammafcn(expon + skew * 0.5j)) - loggammafcn(expon)) - betalnfcn(expon - 0.5, 0.5)
+    return (amplitude / sigma) * exp(logprefactor - expon * log1p(arg * arg) - skew * arctan(arg))
+
+
 def pearson7(x, amplitude=1.0, center=0.0, sigma=1.0, expon=1.0):
     """Return a Pearson7 lineshape.
 

lmfit/models.py

@@ -9,7 +9,7 @@
 from .lineshapes import (breit_wigner, damped_oscillator, dho, doniach,
                          expgaussian, exponential, gaussian, gaussian2d,
                          linear, lognormal, lorentzian, moffat, parabolic,
-                         pearson7, powerlaw, pvoigt, rectangle, sine,
+                         pearson4, pearson7, powerlaw, pvoigt, rectangle, sine,
                          skewed_gaussian, skewed_voigt, split_lorentzian, step,
                          students_t, thermal_distribution, tiny, voigt)
 from .model import Model
@@ -735,6 +735,59 @@ def guess(self, data, x, negative=False, **kwargs):
     guess.__doc__ = COMMON_GUESS_DOC
 
 
+class Pearson4Model(Model):
+    r"""A model based on a Pearson IV distribution.
+
+    The model has five parameters: `amplitude` (:math:`A`), `center`
+    (:math:`\mu`), `sigma` (:math:`\sigma`), `expon` (:math:`m`) and `skew` (:math:`\nu`).
+    In addition, parameters `fwhm`, `height` and `position` are included as
+    constraints to report estimates for the approximate full width at half maximum (20% error),
+    the peak height, and the peak position (the position of the maximal  function value), respectively.
+    The fwhm value has an error of about 20% in the
+    parameter range expon: (0.5, 1000], skew: [-1000, 1000].
+
+    .. math::
+
+        f(x;A,\mu,\sigma,m,\nu)=A \frac{\left|\frac{\Gamma(m+i\tfrac{\nu}{2})}{\Gamma(m)}\right|^2}{\sigma\beta(m-\tfrac{1}{2},\tfrac{1}{2})}\left[1+\frac{(x-\mu)^2}{\sigma^2}\right]^{-m}\exp\left(-\nu \arctan\left(\frac{x-\mu}{\sigma}\right)\right)
+
+    where :math:`\beta` is the beta function (see :scipydoc:`special.beta`).
+    The :meth:`guess` function always gives a starting value of 1.5 for `expon`,
+    and 0 for `skew`.
+
+    For more information, see:
+    https://en.wikipedia.org/wiki/Pearson_distribution#The_Pearson_type_IV_distribution
+
+    """
+
+    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
+                 **kwargs):
+        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
+                       'independent_vars': independent_vars})
+        super().__init__(pearson4, **kwargs)
+        self._set_paramhints_prefix()
+
+    def _set_paramhints_prefix(self):
+        self.set_param_hint('expon', value=1.5, min=0.5 + tiny, max=1000)
+        self.set_param_hint('skew', value=0.0, min=-1000, max=1000)
+        fmt = ("{prefix:s}sigma*sqrt(2**(1/{prefix:s}expon)-1)*pi/arctan2(exp(1)*{prefix:s}expon, {prefix:s}skew)")
+        self.set_param_hint('fwhm', expr=fmt.format(prefix=self.prefix))
+        fmt = ("({prefix:s}amplitude / {prefix:s}sigma) * exp(2 * (real(loggammafcn({prefix:s}expon + {prefix:s}skew * 0.5j)) - loggammafcn({prefix:s}expon)) - betalnfnc({prefix:s}expon-0.5, 0.5) - "
+               "{prefix:s}expon * log1p(square({prefix:s}skew/(2*{prefix:s}expon))) - {prefix:s}skew * arctan(-{prefix:s}skew/(2*{prefix:s}expon)))")
+        self.set_param_hint('height', expr=fmt.format(tiny, prefix=self.prefix))
+        fmt = ("{prefix:s}center-{prefix:s}sigma*{prefix:s}skew/(2*{prefix:s}expon)")
+        self.set_param_hint('position', expr=fmt.format(prefix=self.prefix))
+
+    def guess(self, data, x, negative=False, **kwargs):
+        """Estimate initial model parameter values from data."""
+        pars = guess_from_peak(self, data, x, negative)
+        pars[f'{self.prefix}expon'].set(value=1.5)
+        pars[f'{self.prefix}skew'].set(value=0.0)
+        return update_param_vals(pars, self.prefix, **kwargs)
+
+    __init__.__doc__ = COMMON_INIT_DOC
+    guess.__doc__ = COMMON_GUESS_DOC
+
+
 class Pearson7Model(Model):
     r"""A model based on a Pearson VII distribution.
 
@@ -1531,6 +1584,7 @@ def _parse_params(self):
                 'Voigt': VoigtModel,
                 'PseudoVoigt': PseudoVoigtModel,
                 'Moffat': MoffatModel,
+                'Pearson4': Pearson4Model,
                 'Pearson7': Pearson7Model,
                 'StudentsT': StudentsTModel,
                 'Breit-Wigner': BreitWignerModel,

lmfit/parameter.py

@@ -11,7 +11,7 @@
 from .lineshapes import tiny
 from .printfuncs import params_html_table
 
-SCIPY_FUNCTIONS = {'gamfcn': scipy.special.gamma}
+SCIPY_FUNCTIONS = {'gamfcn': scipy.special.gamma, 'loggammafcn': scipy.special.loggamma, 'betalnfnc': scipy.special.betaln}
 for fnc_name in ('erf', 'erfc', 'wofz'):
     SCIPY_FUNCTIONS[fnc_name] = getattr(scipy.special, fnc_name)
 

tests/test_builtin_models.py

@@ -20,6 +20,8 @@ def check_height_fwhm(x, y, lineshape, model):
     mu = out.params['center'].value
     if lineshape is lineshapes.lognormal:
         cen = np.exp(mu - out.params['sigma']**2)
+    elif lineshape is lineshapes.pearson4:
+        cen = out.params['position']
     else:
         cen = mu
 
@@ -66,6 +68,7 @@ def test_height_fwhm_calculation(peakdata):
     y = peakdata[1]
     check_height_fwhm(x, y, lineshapes.voigt, models.VoigtModel())
     check_height_fwhm(x, y, lineshapes.pvoigt, models.PseudoVoigtModel())
+    check_height_fwhm(x, y, lineshapes.pearson4, models.Pearson4Model())
     check_height_fwhm(x, y, lineshapes.pearson7, models.Pearson7Model())
     check_height_fwhm(x, y, lineshapes.moffat, models.MoffatModel())
     check_height_fwhm(x, y, lineshapes.students_t, models.StudentsTModel())
@@ -112,6 +115,10 @@ def test_height_and_fwhm_expression_evalution_in_builtin_models():
     params = mod.make_params(amplitude=1.0, center=0.0, sigma=0.9, beta=0.0)
     params.update_constraints()
 
+    mod = models.Pearson4Model()
+    params = mod.make_params(amplitude=1.0, center=0.0, sigma=0.9, expon=1.0, skew=5.0)
+    params.update_constraints()
+
     mod = models.Pearson7Model()
     params = mod.make_params(amplitude=1.0, center=0.0, sigma=0.9, expon=1.0)
     params.update_constraints()
@@ -215,6 +222,12 @@ def test_guess_modelparams():
     assert_allclose(pars['l_center'].value, 0.25, rtol=1)
     assert_allclose(pars['l_sigma'].value, 1.3, rtol=1)
 
+    mod = models.Pearson4Model(prefix='g_')
+    pars = mod.guess(y, x=x)
+    assert_allclose(pars['g_amplitude'].value, 3, rtol=2)
+    assert_allclose(pars['g_center'].value, 0.25, rtol=1)
+    assert_allclose(pars['g_sigma'].value, 1.3, rtol=1)
+
     mod = models.SplitLorentzianModel(prefix='s_')
     pars = mod.guess(y, x=x)
     assert_allclose(pars['s_amplitude'].value, 3, rtol=2)