diff --git a/docs/api_reference.rst b/docs/api_reference.rst index 9f18bb80..79f56106 100644 --- a/docs/api_reference.rst +++ b/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: diff --git a/preliz/distributions/continuous.py b/preliz/distributions/continuous.py index 40b84094..617b2491 100644 --- a/preliz/distributions/continuous.py +++ b/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. diff --git a/preliz/distributions/vonmises.py b/preliz/distributions/vonmises.py new file mode 100644 index 00000000..f5bc1a95 --- /dev/null +++ b/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() diff --git a/preliz/internal/optimization.py b/preliz/internal/optimization.py index c3b0728b..d1d586f6 100644 --- a/preliz/internal/optimization.py +++ b/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 diff --git a/preliz/tests/test_scipy.py b/preliz/tests/test_scipy.py index fb1120eb..af4d8222 100644 --- a/preliz/tests/test_scipy.py +++ b/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")