Add faster Weibull

docs/api_reference.rst

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

preliz/distributions/beta.py

@@ -227,7 +227,7 @@ def nb_ppf(q, alpha, beta, lower, upper):
     return ppf_bounds_cont(x_val, q, lower, upper)
 
 
-# @nb.njit
+@nb.njit
 def nb_entropy(alpha, beta):
     psc = alpha + beta
     return (

preliz/distributions/continuous.py

@@ -14,11 +14,12 @@
 from scipy.special import logit, expit  # pylint: disable=no-name-in-module
 
 from ..internal.optimization import optimize_ml, optimize_moments, optimize_moments_rice
-from ..internal.distribution_helper import garcia_approximation, all_not_none, any_not_none
+from ..internal.distribution_helper import all_not_none, any_not_none
 from .distributions import Continuous
 from .beta import Beta  # pylint: disable=unused-import
 from .normal import Normal  # pylint: disable=unused-import
 from .halfnormal import HalfNormal  # pylint: disable=unused-import
+from .weibull import Weibull  # pylint: disable=unused-import
 
 
 eps = np.finfo(float).eps
@@ -2601,75 +2602,3 @@ def _fit_moments(self, mean, sigma):
     def _fit_mle(self, sample, **kwargs):
         mu, _, lam = self.dist.fit(sample, **kwargs)
         self._update(mu * lam, lam)
-
-
-class Weibull(Continuous):
-    r"""
-    Weibull distribution.
-
-    The pdf of this distribution is
-
-    .. math::
-
-       f(x \mid \alpha, \beta) =
-           \frac{\alpha x^{\alpha - 1}
-           \exp(-(\frac{x}{\beta})^{\alpha})}{\beta^\alpha}
-
-    .. plot::
-        :context: close-figs
-
-        import arviz as az
-        from preliz import Weibull
-        az.style.use('arviz-white')
-        alphas = [1., 2, 5.]
-        betas = [1., 1., 2.]
-        for a, b in zip(alphas, betas):
-            Weibull(a, b).plot_pdf(support=(0,5))
-
-    ========  ====================================================
-    Support   :math:`x \in [0, \infty)`
-    Mean      :math:`\beta \Gamma(1 + \frac{1}{\alpha})`
-    Variance  :math:`\beta^2 \Gamma(1 + \frac{2}{\alpha} - \mu^2/\beta^2)`
-    ========  ====================================================
-
-    Parameters
-    ----------
-    alpha : float
-        Shape parameter (alpha > 0).
-    beta : float
-        Scale parameter (beta > 0).
-    """
-
-    def __init__(self, alpha=None, beta=None):
-        super().__init__()
-        self.dist = copy(stats.weibull_min)
-        self.support = (0, np.inf)
-        self._parametrization(alpha, beta)
-
-    def _parametrization(self, alpha=None, beta=None):
-        self.alpha = alpha
-        self.beta = beta
-        self.param_names = ("alpha", "beta")
-        self.params_support = ((eps, np.inf), (eps, np.inf))
-        if all_not_none(alpha, beta):
-            self._update(alpha, beta)
-
-    def _get_frozen(self):
-        frozen = None
-        if all_not_none(self.params):
-            frozen = self.dist(self.alpha, scale=self.beta)
-        return frozen
-
-    def _update(self, alpha, beta):
-        self.alpha = np.float64(alpha)
-        self.beta = np.float64(beta)
-        self.params = (self.alpha, self.beta)
-        self._update_rv_frozen()
-
-    def _fit_moments(self, mean, sigma):
-        alpha, beta = garcia_approximation(mean, sigma)
-        self._update(alpha, beta)
-
-    def _fit_mle(self, sample, **kwargs):
-        alpha, _, beta = self.dist.fit(sample, **kwargs)
-        self._update(alpha, beta)

preliz/distributions/weibull.py

@@ -0,0 +1,170 @@
+# pylint: disable=attribute-defined-outside-init
+# pylint: disable=arguments-differ
+import numba as nb
+import numpy as np
+
+from .distributions import Continuous
+from ..internal.distribution_helper import eps, all_not_none
+from ..internal.optimization import optimize_ml
+from ..internal.special import (
+    garcia_approximation,
+    gamma,
+    mean_and_std,
+    cdf_bounds,
+    ppf_bounds_cont,
+)
+
+
+class Weibull(Continuous):
+    r"""
+    Weibull distribution.
+
+    The pdf of this distribution is
+
+    .. math::
+
+       f(x \mid \alpha, \beta) =
+           \frac{\alpha x^{\alpha - 1}
+           \exp(-(\frac{x}{\beta})^{\alpha})}{\beta^\alpha}
+
+    .. plot::
+        :context: close-figs
+
+        import arviz as az
+        from preliz import Weibull
+        az.style.use('arviz-white')
+        alphas = [1., 2, 5.]
+        betas = [1., 1., 2.]
+        for a, b in zip(alphas, betas):
+            Weibull(a, b).plot_pdf(support=(0,5))
+
+    ========  ====================================================
+    Support   :math:`x \in [0, \infty)`
+    Mean      :math:`\beta \Gamma(1 + \frac{1}{\alpha})`
+    Variance  :math:`\beta^2 \Gamma(1 + \frac{2}{\alpha} - \mu^2/\beta^2)`
+    ========  ====================================================
+
+    Parameters
+    ----------
+    alpha : float
+        Shape parameter (alpha > 0).
+    beta : float
+        Scale parameter (beta > 0).
+    """
+
+    def __init__(self, alpha=None, beta=None):
+        super().__init__()
+        self.support = (0, np.inf)
+        self._parametrization(alpha, beta)
+
+    def _parametrization(self, alpha=None, beta=None):
+        self.alpha = alpha
+        self.beta = beta
+        self.param_names = ("alpha", "beta")
+        self.params_support = ((eps, np.inf), (eps, np.inf))
+        if all_not_none(alpha, beta):
+            self._update(alpha, beta)
+
+    def _update(self, alpha, beta):
+        self.alpha = np.float64(alpha)
+        self.beta = np.float64(beta)
+        self.params = (self.alpha, self.beta)
+
+        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(nb_logpdf(x, self.alpha, self.beta))
+
+    def cdf(self, x):
+        """
+        Compute the cumulative distribution function (CDF) at a given point x.
+        """
+        x = np.asarray(x)
+        return nb_cdf(x, self.alpha, self.beta, self.support[0], self.support[1])
+
+    def ppf(self, q):
+        """
+        Compute the percent point function (PPF) at a given probability q.
+        """
+        q = np.asarray(q)
+        return nb_ppf(q, self.alpha, self.beta, self.support[0], self.support[1])
+
+    def logpdf(self, x):
+        """
+        Compute the log probability density function (log PDF) at a given point x.
+        """
+        return nb_logpdf(x, self.alpha, self.beta)
+
+    def entropy(self):
+        return nb_entropy(self.alpha, self.beta)
+
+    def mean(self):
+        return self.beta * gamma(1 + 1 / self.alpha)
+
+    def median(self):
+        return self.beta * np.log(2) ** (1 / self.alpha)
+
+    def var(self):
+        return self.beta**2 * gamma(1 + 2 / self.alpha) - self.mean() ** 2
+
+    def std(self):
+        return self.var() ** 0.5
+
+    def skewness(self):
+        mu = self.mean()
+        sigma = self.std()
+        m_s = mu / sigma
+        return gamma(1 + 3 / self.alpha) * (self.beta / sigma) ** 3 - 3 * m_s - m_s**3
+
+    def kurtosis(self):
+        mu = self.mean()
+        sigma = self.std()
+        skew = self.skewness()
+        m_s = mu / sigma
+        return (
+            (self.beta / sigma) ** 4 * gamma(1 + 4 / self.alpha)
+            - 4 * skew * m_s
+            - 6 * m_s**2
+            - m_s**4
+            - 3
+        )
+
+    def rvs(self, size=1, random_state=None):
+        random_state = np.random.default_rng(random_state)
+        return random_state.weibull(self.alpha, size) * self.beta
+
+    def _fit_moments(self, mean, sigma):
+        alpha, beta = garcia_approximation(mean, sigma)
+        self._update(alpha, beta)
+
+    def _fit_mle(self, sample):
+        mean, std = mean_and_std(sample)
+        self._fit_moments(mean, std)
+        optimize_ml(self, sample)
+
+
+@nb.njit
+def nb_cdf(x, alpha, beta, lower, upper):
+    prob = 1 - np.exp(-((x / beta) ** alpha))
+    return cdf_bounds(prob, x, lower, upper)
+
+
+@nb.njit
+def nb_ppf(q, alpha, beta, lower, upper):
+    x_val = beta * (-np.log(1 - q)) ** (1 / alpha)
+    return ppf_bounds_cont(x_val, q, lower, upper)
+
+
+@nb.njit
+def nb_entropy(alpha, beta):
+    return np.euler_gamma * (1 - 1 / alpha) + np.log(beta / alpha) + 1
+
+
+@nb.njit
+def nb_logpdf(x, alpha, beta):
+    x_b = x / beta
+    return np.log(alpha / beta) + (alpha - 1) * np.log(x_b) - x_b**alpha

preliz/internal/distribution_helper.py

@@ -1,7 +1,5 @@
 import re
 import numpy as np
-from scipy.special import gamma
-
 
 eps = np.finfo(float).eps
 
@@ -43,27 +41,6 @@ def hdi_from_pdf(dist, mass=0.94):
     return x_vals[np.sort(indices)[[0, -1]]]
 
 
-def garcia_approximation(mean, sigma):
-    """
-    Approximate method of moments for Weibull distribution, provides good results for values of
-    alpha larger than 0.83.
-
-    Oscar Garcia. Simplified method-of-moments estimation for the Weibull distribution. 1981.
-    New Zealand Journal of Forestry Science 11:304-306
-    https://www.scionresearch.com/__data/assets/pdf_file/0010/36874/NZJFS1131981GARCIA304-306.pdf
-    """
-    ks = [-0.221016417, 0.010060668, 0.117358987, -0.050999126]  # pylint: disable=invalid-name
-    z = sigma / mean  # pylint: disable=invalid-name
-
-    poly = 0
-    for idx, k in enumerate(ks):
-        poly += k * z**idx
-
-    alpha = 1 / (z * (1 + (1 - z) ** 2 * poly))
-    beta = 1 / (gamma(1 + 1 / alpha) / (mean))
-    return alpha, beta
-
-
 def all_not_none(*args):
     for arg in args:
         if arg is None: