Add faster VonMises

docs/api_reference.rst

@@ -42,6 +42,9 @@ This reference provides detailed documentation for user functions in the current
 .. automodule:: preliz.distributions.normal
    :members:
 
+.. automodule:: preliz.distributions.vonmises
+   :members:
+
 .. automodule:: preliz.distributions.weibull
    :members:
 

preliz/distributions/continuous.py

@@ -2,6 +2,7 @@
 # pylint: disable=too-many-instance-attributes
 # pylint: disable=invalid-name
 # pylint: disable=attribute-defined-outside-init
+# pylint: disable=unused-import
 """
 Continuous probability distributions.
 """
@@ -16,13 +17,14 @@
 from ..internal.optimization import optimize_ml, optimize_moments, optimize_moments_rice
 from ..internal.distribution_helper import all_not_none, any_not_none
 from .distributions import Continuous
-from .asymmetric_laplace import AsymmetricLaplace  # pylint: disable=unused-import
-from .beta import Beta  # pylint: disable=unused-import
-from .exponential import Exponential  # pylint: disable=unused-import
-from .normal import Normal  # pylint: disable=unused-import
-from .halfnormal import HalfNormal  # pylint: disable=unused-import
-from .laplace import Laplace  # pylint: disable=unused-import
-from .weibull import Weibull  # pylint: disable=unused-import
+from .asymmetric_laplace import AsymmetricLaplace
+from .beta import Beta
+from .exponential import Exponential
+from .normal import Normal
+from .halfnormal import HalfNormal
+from .laplace import Laplace
+from .vonmises import VonMises
+from .weibull import Weibull
 
 
 eps = np.finfo(float).eps
@@ -2130,80 +2132,6 @@ def _fit_mle(self, sample, **kwargs):
         self._update(lower, upper)
 
 
-class VonMises(Continuous):
-    r"""
-    Univariate VonMises distribution.
-
-    The pdf of this distribution is
-
-    .. math::
-
-        f(x \mid \mu, \kappa) =
-            \frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)}
-
-    where :math:`I_0` is the modified Bessel function of order 0.
-
-    .. plot::
-        :context: close-figs
-
-        import arviz as az
-        from preliz import VonMises
-        az.style.use('arviz-doc')
-        mus = [0., 0., 0.,  -2.5]
-        kappas = [.01, 0.5, 4., 2.]
-        for mu, kappa in zip(mus, kappas):
-            VonMises(mu, kappa).plot_pdf(support=(-np.pi,np.pi))
-
-    ========  ==========================================
-    Support   :math:`x \in [-\pi, \pi]`
-    Mean      :math:`\mu`
-    Variance  :math:`1-\frac{I_1(\kappa)}{I_0(\kappa)}`
-    ========  ==========================================
-
-    Parameters
-    ----------
-    mu : float
-        Mean.
-    kappa : float
-        Concentration (:math:`\frac{1}{\kappa}` is analogous to :math:`\sigma^2`).
-    """
-
-    def __init__(self, mu=None, kappa=None):
-        super().__init__()
-        self.dist = copy(stats.vonmises)
-        self._parametrization(mu, kappa)
-
-    def _parametrization(self, mu=None, kappa=None):
-        self.mu = mu
-        self.kappa = kappa
-        self.param_names = ("mu", "kappa")
-        self.params_support = ((-np.pi, np.pi), (eps, np.inf))
-        self.support = (-np.pi, np.pi)
-        if all_not_none(mu, kappa):
-            self._update(mu, kappa)
-
-    def _get_frozen(self):
-        frozen = None
-        if all_not_none(self.params):
-            frozen = self.dist(kappa=self.kappa, loc=self.mu)
-        return frozen
-
-    def _update(self, mu, kappa):
-        self.mu = np.float64(mu)
-        self.kappa = np.float64(kappa)
-        self.params = (self.mu, self.kappa)
-        self._update_rv_frozen()
-
-    def _fit_moments(self, mean, sigma):
-        mu = mean
-        kappa = 1 / sigma**2
-        self._update(mu, kappa)
-
-    def _fit_mle(self, sample, **kwargs):
-        kappa, mu, _ = self.dist.fit(sample, **kwargs)
-        self._update(mu, kappa)
-
-
 class Wald(Continuous):
     r"""
     Wald distribution.

preliz/distributions/vonmises.py

@@ -0,0 +1,187 @@
+# pylint: disable=attribute-defined-outside-init
+# pylint: disable=arguments-differ
+import numpy as np
+from scipy.special import i0e, i1e  # pylint: disable=no-name-in-module
+from scipy.integrate import quad
+from scipy.optimize import bisect
+from scipy.stats import circmean
+from .distributions import Continuous
+from ..internal.distribution_helper import eps, all_not_none
+from ..internal.optimization import find_kappa, optimize_moments
+
+
+class VonMises(Continuous):
+    r"""
+    Univariate VonMises distribution.
+
+    The pdf of this distribution is
+
+    .. math::
+
+        f(x \mid \mu, \kappa) =
+            \frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)}
+
+    where :math:`I_0` is the modified Bessel function of order 0.
+
+    .. plot::
+        :context: close-figs
+
+        import arviz as az
+        from preliz import VonMises
+        az.style.use('arviz-doc')
+        mus = [0., 0., 0.,  -2.5]
+        kappas = [.01, 0.5, 4., 2.]
+        for mu, kappa in zip(mus, kappas):
+            VonMises(mu, kappa).plot_pdf(support=(-np.pi,np.pi))
+
+    ========  ==========================================
+    Support   :math:`x \in [-\pi, \pi]`
+    Mean      :math:`\mu`
+    Variance  :math:`1-\frac{I_1(\kappa)}{I_0(\kappa)}`
+    ========  ==========================================
+
+    Parameters
+    ----------
+    mu : float
+        Mean.
+    kappa : float
+        Concentration (:math:`\frac{1}{\kappa}` is analogous to :math:`\kappa^2`).
+    """
+
+    def __init__(self, mu=None, kappa=None):
+        super().__init__()
+        self._parametrization(mu, kappa)
+
+    def _parametrization(self, mu=None, kappa=None):
+        self.mu = mu
+        self.kappa = kappa
+        self.param_names = ("mu", "kappa")
+        self.params_support = ((-np.pi, np.pi), (eps, np.inf))
+        self.support = (-np.pi, np.pi)
+        if all_not_none(mu, kappa):
+            self._update(mu, kappa)
+
+    def _update(self, mu, kappa):
+        self.mu = np.float64(mu)
+        self.kappa = np.float64(kappa)
+        self.params = (self.mu, self.kappa)
+        self.is_frozen = True
+
+    def pdf(self, x):
+        """
+        Compute the probability density function (PDF) at a given point x.
+        """
+        x = np.asarray(x)
+        return np.exp(self.logpdf(x))
+
+    def cdf(self, x):
+        """
+        Compute the cumulative distribution function (CDF) at a given point x.
+        """
+        return nb_cdf(x, self.pdf)
+
+    def ppf(self, q):
+        """
+        Compute the percent point function (PPF) at a given probability q.
+        """
+        return nb_ppf(q, self.pdf)
+
+    def logpdf(self, x):
+        """
+        Compute the log probability density function (log PDF) at a given point x.
+        """
+        return nb_logpdf(x, self.mu, self.kappa)
+
+    def _neg_logpdf(self, x):
+        """
+        Compute the neg log_pdf sum for the array x.
+        """
+        return nb_neg_logpdf(x, self.mu, self.kappa)
+
+    def entropy(self):
+        return nb_entropy(self.kappa, self.var())
+
+    def mean(self):
+        return self.mu
+
+    def median(self):
+        return self.mu
+
+    def var(self):
+        return 1 - i1e(self.kappa) / i0e(self.kappa)
+
+    def std(self):
+        return self.var() ** 0.5
+
+    def skewness(self):
+        return 0
+
+    def kurtosis(self):
+        return 0
+
+    def rvs(self, size=1, random_state=None):
+        random_state = np.random.default_rng(random_state)
+        return random_state.vonmises(self.mu, self.kappa, size)
+
+    def _fit_moments(self, mean, sigma):
+        params = mean, 1 / sigma**1.8
+        optimize_moments(self, mean, sigma, params)
+
+    def _fit_mle(self, sample):
+        data = np.mod(sample, 2 * np.pi)
+        mu = circmean(data)
+        kappa = find_kappa(data, mu)
+        mu = np.mod(mu + np.pi, 2 * np.pi) - np.pi
+        self._update(mu, kappa)
+
+
+def nb_cdf(x, pdf):
+    if isinstance(x, (int, float)):
+        x = [x]
+        scalar_input = True
+    else:
+        scalar_input = False
+
+    cdf_values = np.array([quad(pdf, -np.pi, xi)[0] if xi <= np.pi else 1 for xi in x])
+
+    return cdf_values[0] if scalar_input else cdf_values
+
+
+def nb_ppf(q, pdf):
+    def root_func(x, q):
+        return nb_cdf(x, pdf) - q
+
+    if isinstance(q, (int, float)):
+        q = [q]
+        scalar_input = True
+    else:
+        scalar_input = False
+
+    ppf_values = []
+    for q_i in q:
+        if q_i < 0:
+            val = np.nan
+        elif q_i > 1:
+            val = np.nan
+        elif q_i == 0:
+            val = -np.inf
+        elif q_i == 1:
+            val = np.inf
+        else:
+            val = bisect(root_func, -np.pi, np.pi, args=(q_i,))
+
+        ppf_values.append(val)
+
+    return ppf_values[0] if scalar_input else np.array(ppf_values)
+
+
+def nb_entropy(kappa, var):
+    return np.log(2 * np.pi * i0e(kappa)) + kappa * var
+
+
+def nb_logpdf(x, mu, kappa):
+    return kappa * (np.cos(x - mu) - 1) - np.log(2 * np.pi) - np.log(i0e(kappa))
+
+
+def nb_neg_logpdf(x, mu, kappa):
+    return -(nb_logpdf(x, mu, kappa)).sum()

preliz/internal/optimization.py

@@ -6,8 +6,8 @@
 from copy import copy
 
 import numpy as np
-from scipy.optimize import minimize, least_squares
-from scipy.special import i0, i1  # pylint: disable=no-name-in-module
+from scipy.optimize import minimize, least_squares, root_scalar
+from scipy.special import i0, i1, i0e, i1e  # pylint: disable=no-name-in-module
 from .distribution_helper import init_vals as default_vals
 
 
@@ -397,3 +397,19 @@ def get_fixed_params(distribution):
         else:
             fixed.append(value)
     return none_idx, fixed
+
+
+def find_kappa(data, mu):
+    ere = np.mean(np.cos(mu - data))
+
+    if ere > 0:
+
+        def solve_for_kappa(kappa):
+            return i1e(kappa) / i0e(kappa) - ere
+
+        root_res = root_scalar(
+            solve_for_kappa, method="brentq", bracket=(np.finfo(float).tiny, 1e16)
+        )
+        return root_res.root
+    else:
+        return np.finfo(float).tiny

preliz/tests/test_scipy.py

@@ -11,6 +11,7 @@
     HalfNormal,
     Laplace,
     Normal,
+    VonMises,
     Weibull,
     Bernoulli,
     Binomial,
@@ -36,6 +37,7 @@
         (HalfNormal, stats.halfnorm, {"sigma": 2}, {"scale": 2}),
         (Laplace, stats.laplace, {"mu": 2.5, "b": 4}, {"loc": 2.5, "scale": 4}),
         (Normal, stats.norm, {"mu": 0, "sigma": 2}, {"loc": 0, "scale": 2}),
+        (VonMises, stats.vonmises, {"mu": 0, "kappa": 10}, {"loc": 0, "kappa": 10}),
         (
             Weibull,
             stats.weibull_min,
@@ -75,12 +77,15 @@ def test_match_scipy(p_dist, sp_dist, p_params, sp_params):
     preliz_dist = p_dist(**p_params)
     scipy_dist = sp_dist(**sp_params)
 
-    actual = preliz_dist.entropy()
-    expected = scipy_dist.entropy()
-    if preliz_dist.kind == "discrete":
-        assert_almost_equal(actual, expected, decimal=1)
-    else:
-        assert_almost_equal(actual, expected, decimal=4)
+    if preliz_dist.__class__.__name__ != "VonMises":
+        # for the VonMises we used the differential entropy definition.
+        # SciPy uses a different one
+        actual = preliz_dist.entropy()
+        expected = scipy_dist.entropy()
+        if preliz_dist.kind == "discrete":
+            assert_almost_equal(actual, expected, decimal=1)
+        else:
+            assert_almost_equal(actual, expected, decimal=4)
 
     rng = np.random.default_rng(1)
     actual_rvs = preliz_dist.rvs(20, random_state=rng)
@@ -123,9 +128,18 @@ def test_match_scipy(p_dist, sp_dist, p_params, sp_params):
         "ZeroInflatedBinomial",
         "ZeroInflatedNegativeBinomial",
         "ZeroInflatedPoisson",
+        "VonMises",
     ]:
         actual_moments = preliz_dist.moments("mvsk")
         expected_moments = scipy_dist.stats("mvsk")
+    elif preliz_dist.__class__.__name__ == "VonMises":
+        # We use the circular variance definition for the variance
+        assert_almost_equal(preliz_dist.var(), stats.circvar(preliz_dist.rvs(1000)), decimal=1)
+        # And we adopt the convention of setting the skewness and kurtosis to 0 in
+        # analogy with the Normal distribution
+        actual_moments = preliz_dist.moments("m")
+        expected_moments = scipy_dist.stats("m")
+
     else:
         actual_moments = preliz_dist.moments("mv")
         expected_moments = scipy_dist.stats("mv")