Add mode method to more distribution classes

preliz/distributions/continuous_multivariate.py

@@ -307,6 +307,7 @@ def plot(**args):
                     interval,
                     levels,
                     "full",
+                    legend,
                     figsize,
                     None,
                     xy_lim,
@@ -320,13 +321,23 @@ def plot(**args):
                     interval,
                     levels,
                     None,
+                    legend,
                     figsize,
                     None,
                     xy_lim,
                 )
 
         return interactive(plot, **plot_widgets)
 
+    def mode(self):
+        alpha_sum = np.sum(self.alpha)
+        K = len(self.alpha)
+        return np.where(
+            np.all(self.alpha > 1),
+            (self.alpha - 1) / (alpha_sum - K),
+            np.nan
+        )
+
 
 class MvNormal(Continuous):
     r"""
@@ -669,3 +680,10 @@ def plot(**args):
                 )
 
         return interactive(plot, **plot_widgets)
+
+    def mode(self):
+        """
+        Calculate the mode of the Multivariate Normal distribution.
+        For Multivariate Normal, the mode equals the mean.
+        """
+        return self.mu

preliz/distributions/discrete_weibull.py

@@ -5,7 +5,7 @@
 from preliz.internal.distribution_helper import all_not_none, eps, num_kurtosis, num_skewness
 from preliz.internal.optimization import optimize_ml, optimize_moments
 from preliz.internal.special import cdf_bounds, ppf_bounds_disc
-
+from scipy.optimize import minimize_scalar
 
 class DiscreteWeibull(Discrete):
     R"""
@@ -107,6 +107,14 @@ def skewness(self):
     def kurtosis(self):
         return num_kurtosis(self)
 
+    def mode(self):
+
+        def negative_pdf(x):
+            return -self.pdf(x)
+
+        result = minimize_scalar(negative_pdf, bounds=(0, 100), method='bounded')
+        return np.floor(result.x)
+
     def rvs(self, size=None, random_state=None):
         random_state = np.random.default_rng(random_state)
         return self.ppf(random_state.uniform(size=size))

preliz/distributions/exgaussian.py

@@ -1,6 +1,7 @@
 import numba as nb
 import numpy as np
 from scipy.stats import skew
+from scipy.special import erfcinv
 
 from preliz.distributions.distributions import Continuous
 from preliz.internal.distribution_helper import all_not_none, eps
@@ -101,6 +102,14 @@ def entropy(self):
         logpdf = self.logpdf(x_values)
         return -np.trapz(np.exp(logpdf) * logpdf, x_values)
 
+    def mode(self):
+        tau = self.nu
+        mu = self.mu
+        sigma = self.sigma
+
+        t = np.abs(tau) * np.sqrt(2/np.pi)
+        return mu - np.sign(tau) * np.sqrt(2*sigma) * erfcinv(t) + sigma**2/tau
+
     def mean(self):
         return self.mu + self.nu
 

preliz/distributions/hypergeometric.py

@@ -143,6 +143,9 @@ def kurtosis(self):
             )
         )
 
+    def mode(self):
+        return np.floor((self.n + 1) * (self.k + 1) / (self.N + 2))
+
     def rvs(self, size=None, random_state=None):
         random_state = np.random.default_rng(random_state)
         return random_state.hypergeometric(self.k, self.N - self.k, self.n, size=size)

preliz/distributions/logitnormal.py

@@ -1,6 +1,6 @@
 import numba as nb
 import numpy as np
-
+from scipy.optimize import root_scalar
 from preliz.distributions.distributions import Continuous
 from preliz.internal.distribution_helper import all_not_none, eps, from_precision, to_precision
 from preliz.internal.special import (
@@ -149,6 +149,37 @@ def kurtosis(self):
         pdf = self.pdf(x_values)
         return np.trapz(((x_values - mean) / std) ** 4 * pdf, x_values) - 3
 
+    def mode(self):
+        def mode_equation(x):
+            # The equation is: logit(x) = σ²(2x-1) + μ
+            # We want to find the root of: logit(x) - σ²(2x-1) - μ = 0
+            return logit(x) - (self.sigma**2 * (2*x - 1)) - self.mu     
+
+        #Left side
+        try:
+            sol1 = root_scalar(mode_equation, bracket=(eps, 0.5-eps)).root
+        except ValueError:
+            sol1 = None
+
+        #Right side
+        try:
+            sol2 = root_scalar(mode_equation, bracket=(0.5+eps, 1-eps)).root
+        except ValueError:
+            sol2 = None
+
+        if sol1 is None and sol2 is None:
+            # If no solutions found, return the median as an approximation
+            return self.median()
+        elif sol1 is None:
+            return sol2
+        elif sol2 is None:
+            return sol1
+        else:
+            # Return the solution with higher density
+            if self.pdf(sol1) >= self.pdf(sol2):
+                return sol1
+            return sol2
+
     def rvs(self, size=None, random_state=None):
         random_state = np.random.default_rng(random_state)
         return expit(random_state.normal(self.mu, self.sigma, size))

preliz/distributions/skewnormal.py

@@ -150,6 +150,21 @@ def kurtosis(self):
             * ((delta * np.sqrt(2 / np.pi)) ** 4 / (1 - 2 * (delta**2) / np.pi) ** 2)
         )
 
+    def mode(self):
+        alpha = self.alpha
+        delta = alpha / np.sqrt(1 + alpha**2)
+
+        # Calculate mo(alpha)
+        sqrt_2_pi = np.sqrt(2/np.pi)
+        term1 = sqrt_2_pi * delta
+        term2 = (1 - np.pi/4) * (sqrt_2_pi * delta)**3 / (1 - 2/np.pi * delta**2)
+        term3 = np.sign(alpha)/2 * np.exp(-2*np.pi/abs(alpha)) if alpha != 0 else 0
+
+        mo_alpha = term1 - term2 - term3
+
+        # Final mode calculation
+        return self.mu + self.sigma * mo_alpha
+
     def rvs(self, size=None, random_state=None):
         random_state = np.random.default_rng(random_state)
         u_0 = random_state.normal(size=size)