Addition of Johnson SU pdf (#92)

* Introduction of JohnsonSU pdf.

* addition of johnsonSU plots for doc

* remove johnsonSU from pdf.pxd

* fix missing comma in __init__.py

* finish implementation of normalized johnson pdf and add tests.

* fix error in __init__

* replace .info to .version in probfit/__init__.py (bad previous merging)

* dd in the docstring a reference for the JohnsonSU PDF

doc/api.rst

@@ -221,6 +221,13 @@ poly3
 .. plot:: pyplots/pdf/poly3.py
     :class: lightbox
 
+JohnsonSU
+^^^^^^^^^
+.. autofunction:: johnsonSU
+
+.. plot:: pyplots/pdf/johnsonSU.py
+    :class: lightbox
+
 Polynomial
 ^^^^^^^^^^
 .. autoclass:: Polynomial

doc/index.rst

@@ -104,6 +104,7 @@ and send us pull request.
     cauchy
     rtv_breitwigner
     doublegaussian
+    johnsonSU
     argus
     linear
     poly2

doc/pyplots/pdf/johnsonSU.py

@@ -0,0 +1,21 @@
+from probfit.pdf import johnsonSU
+from probfit.plotting import draw_pdf
+from collections import OrderedDict
+import matplotlib.pyplot as plt
+
+bound = (-10, 10)
+
+arg = OrderedDict(mean=2, sigma=1, nu=-4, tau=0.5)
+draw_pdf(johnsonSU, arg=arg, bound=bound, label=str(arg), density=False,
+         bins=200)
+
+arg = OrderedDict(mean=-3, sigma=2, nu=+4, tau=0.1)
+draw_pdf(johnsonSU, arg=arg, bound=bound, label=str(arg), density=False,
+         bins=200)
+
+arg = OrderedDict(mean=0, sigma=3, nu=+2, tau=0.9)
+draw_pdf(johnsonSU, arg=arg, bound=bound, label=str(arg), density=False,
+         bins=200)
+
+plt.grid(True)
+plt.legend().get_frame().set_alpha(0.5)

probfit/__init__.py

@@ -5,10 +5,11 @@
 * Docs: http://probfit.readthedocs.io
 """
 
-from .costfunc import UnbinnedLH, BinnedLH, Chi2Regression, BinnedChi2, SimultaneousFit
+from .costfunc import UnbinnedLH, BinnedLH, Chi2Regression, BinnedChi2,\
+ SimultaneousFit
 from .pdf import doublegaussian, ugaussian, gaussian, crystalball, \
     argus, cruijff, linear, poly2, poly3, novosibirsk, \
-    Polynomial, HistogramPdf, cauchy, rtv_breitwigner
+    Polynomial, HistogramPdf, cauchy, rtv_breitwigner, johnsonSU
 from .toy import gen_toy, gen_toyn
 from .util import *
 from .oneshot import *
@@ -17,7 +18,9 @@
 from .funcutil import *
 from .decorator import *
 from ._libstat import integrate1d
-from .functor import Normalized, Extended, Convolve, AddPdf, AddPdfNorm, BlindFunc
+
+from .functor import Normalized, Extended, Convolve, AddPdf, AddPdfNorm, \
+ BlindFunc
 from .version import __version__
 
 __all__ = [
@@ -40,6 +43,7 @@
     'crystalball',
     'describe',
     'doublegaussian',
+    'johnsonSU',
     'draw_compare',
     'draw_compare_hist',
     'draw_pdf',

probfit/pdf.pyx

@@ -1,7 +1,7 @@
 #cython: embedsignature=True
 cimport cython
 
-from libc.math cimport exp, pow, fabs, log, sqrt, sinh, tgamma, abs, fabs
+from libc.math cimport exp, pow, fabs, log, sqrt, sinh, tgamma, abs, fabs, cosh, atan2, asinh, erf
 cdef double pi = 3.14159265358979323846264338327
 import numpy as np
 cimport numpy as np
@@ -114,6 +114,100 @@ cdef class HistogramPdf:
             return 0.0
 
 
+cdef class _JohnsonSU:
+    """
+    Normalized JohnsonSU [1]_.
+
+    .. math::
+        f(x; \\mu, \\sigma, \\nu, \\tau) = \\frac{1}{\\lambda \\sqrt{2\\pi}}
+        \\frac{1}{\\sqrt{1 + \\left( \\frac{x - \\xi}{\\lambda} \\right)}}
+        e^{-\\frac{1}{2} \\left( -\\nu + \\frac{1}{\\tau}  \\right)
+        \\sinh^{-1} \\left( \\frac{x - \\xi}{\\lambda} \\right)}
+
+    where
+
+    .. math::
+        \\lambda = \\sigma \\times \\left(\\frac{1}{2}( \\exp(\\tau^{2}) - 1)
+        \\left(\\exp(\\tau^{2}) \\cosh\\left(-\\nu \\tau \\right) + 1\\right)
+        \\right)^{-\\frac{1}{2}} \\\\
+
+    and
+
+    .. math::
+        \\xi = \\mu + \\lambda \\exp\\left(\\frac{\\tau^{2}}{2}\\right)\\sinh
+        \\left( \\nu \\tau \\right)
+
+    References
+    ----------
+
+    .. [1] https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution
+
+    """
+
+    cdef public object func_code
+    cdef public object func_defaults
+
+    def __init__(self, xname='x'):
+
+        varnames = [xname, "mean", "sigma", "nu", "tau"]
+        self.func_code = MinimalFuncCode(varnames)
+        self.func_defaults = None
+
+    def integrate(self, tuple bound, int nint_subdiv=0, *arg):
+
+        cdef double a, b
+        a, b = bound
+
+        cdef double mean = arg[0]
+        cdef double sigma = arg[1]
+        cdef double nu = arg[2]
+        cdef double tau = arg[3]
+
+        cdef double w = exp(tau * tau)
+        cdef double omega = - nu * tau
+        cdef double c = 0.5 * (w-1) * (w * cosh(2 * omega) + 1)
+        c = pow(c, -0.5)
+
+        cdef double _lambda = sigma * c
+        cdef double xi = mean + _lambda * sqrt(w) * sinh(omega)
+        cdef double zmax = (b - xi) / _lambda
+        cdef double zmin = (a - xi) / _lambda
+        cdef double rmax = -nu + asinh(zmax) / tau
+        cdef double rmin = -nu + asinh(zmin) / tau
+
+        cdef double ret = 0.
+
+        ret = 0.5 * (erf(zmax / (sqrt(2))) - erf(zmin / (sqrt(2))))
+
+        return ret
+
+    def __call__(self, *arg):
+
+        cdef double x = arg[0]
+        cdef double mean = arg[1]
+        cdef double sigma = arg[2]
+        cdef double nu = arg[3]
+        cdef double tau = arg[4]
+
+        cdef double w = exp(tau * tau)
+        cdef double omega = - nu * tau
+
+        cdef double c = 0.5 * (w-1) * (w * cosh(2 * omega) + 1)
+        c = pow(c, -0.5)
+
+        cdef double _lambda = sigma * c
+        cdef double xi = mean + _lambda * sqrt(w) * sinh(omega)
+        cdef double z = (x - xi) / _lambda
+        cdef double r = -nu + asinh(z) / tau
+
+        cdef double ret = 1. / (tau * _lambda * sqrt(2 * pi))
+        ret *= 1. / sqrt(z * z + 1)
+        ret *= exp(-0.5 * r * r)
+
+        return ret
+
+johnsonSU = _JohnsonSU()
+
 cpdef double doublegaussian(double x, double mean, double sigma_L, double sigma_R):
     """
     Unnormalized double gaussian

tests/test_func.py

@@ -142,6 +142,18 @@ def test_cauchy():
     assert_allclose(pdf.cauchy(1, 1, 2.), 0.15915494309189535)
     assert_allclose(pdf.cauchy(1, 2, 4.), 0.07489644380795074)
 
+def test_johnsonSU():
+    assert describe(pdf.johnsonSU), ['x', "mean", "sigma", "nu", "tau"]
+    assert_allclose(pdf.johnsonSU(1., 1., 1., 1., 1.), 0.5212726124342)
+    assert_allclose(pdf.johnsonSU(1., 2., 1., 1., 1.), 0.1100533373219)
+    assert_allclose(pdf.johnsonSU(1., 2., 2., 1., 1.), 0.4758433826682)
+
+    j = pdf.johnsonSU
+    assert (hasattr(j, 'integrate'))
+    integral = j.integrate((-100, 100), 0, 1., 1., 1., 1.)
+    assert_allclose(integral, 1.0)
+    integral = j.integrate((0, 2), 0, 1., 1., 1., 1.)
+    assert_allclose(integral, 0.8786191859)
 
 def test_HistogramPdf():
     be = np.array([0, 1, 3, 4], dtype=float)