如何在scipy.stats中获取拟合分布的众数？

Great question! I’ve faced this exact scenario before when working with fitted distributions in SciPy—let’s walk through the two main approaches to find the mode of a fitted distribution.

Core Idea #

First, remember: the mode of a distribution is the value where its probability density function (PDF, for continuous distributions) or probability mass function (PMF, for discrete distributions) reaches its maximum. SciPy doesn’t have a built-in mode() method for fitted distributions, but we can calculate it either using analytical formulas (when available) or numerical optimization.

1. Use Analytical Formulas (For Distributions With Known Mode) #

Many common distributions have a closed-form formula for their mode. You just need to map the fitted parameters from SciPy to the formula. Here are some practical examples:

Example 1: Normal Distribution #

For a normal distribution, the mode equals the mean (and median). After fitting, the loc parameter corresponds directly to this value:

import scipy.stats as stats
import numpy as np

# Generate sample data and fit the distribution
data = np.random.normal(loc=5, scale=2, size=1000)
loc, scale = stats.norm.fit(data)

# Normal distribution mode = loc (mean)
mode = loc
print(f"Normal distribution mode: {round(mode, 2)}")

Example 2: Gamma Distribution #

The mode of a gamma distribution is loc + (shape - 1) * scale when shape > 1. If shape ≤ 1, the PDF is strictly decreasing, so the mode sits at the left boundary (loc):

data = np.random.gamma(shape=3, scale=2, size=1000)
shape, loc, scale = stats.gamma.fit(data)

if shape > 1:
    mode = loc + (shape - 1) * scale
else:
    mode = loc  # PDF starts at loc and decreases from there

print(f"Gamma distribution mode: {round(mode, 2)}")

Example 3: Poisson Distribution (Discrete) #

For a Poisson distribution, the mode is floor(λ) if λ is not an integer. If λ is an integer, both λ and λ-1 are valid modes:

data = np.random.poisson(lam=4.2, size=1000)
lam = stats.poisson.fit(data)[0]

# Check for integer lambda
if lam.is_integer():
    mode = [int(lam - 1), int(lam)]
else:
    mode = int(np.floor(lam))

print(f"Poisson distribution mode: {mode}")

Pro tip: You can find analytical mode formulas for any SciPy distribution in its official docs—look for the "Notes" or "Statistics" section.

2. Numerical Optimization (For Distributions Without Analytical Mode) #

If you don’t remember the formula or the distribution has no closed-form mode, you can numerically find the maximum of the PDF/PMF. For continuous distributions, we minimize the negative PDF (since optimization libraries typically focus on finding minima):

Example: Beta Distribution #

from scipy.optimize import minimize

data = np.random.beta(a=2.5, b=3.5, size=1000)
a, b, loc, scale = stats.beta.fit(data)

# Define the negative PDF function (minimizing this finds the PDF's maximum)
def neg_pdf(x):
    # Ensure x stays within the distribution's valid range (Beta is [loc, loc+scale])
    if not (loc <= x <= loc + scale):
        return np.inf
    return -stats.beta.pdf(x, a, b, loc=loc, scale=scale)

# Start with an initial guess (use the distribution's mean as a starting point)
initial_guess = stats.beta.mean(a, b, loc=loc, scale=scale)

# Run the minimization, with bounds to stay within valid values
result = minimize(neg_pdf, x0=initial_guess, bounds=[(loc, loc+scale)])

if result.success:
    mode = result.x[0]
    print(f"Beta distribution mode (numerical): {round(mode, 2)}")
else:
    print("Failed to compute mode numerically—check your initial guess or bounds.")

For discrete distributions, compute the PMF for a range of plausible values and pick the one(s) with the highest probability:

# Example: Binomial Distribution
data = np.random.binomial(n=10, p=0.6, size=1000)
n, p = stats.binom.fit(data)

# Test all possible values (0 to n)
possible_values = np.arange(0, n+1)
pmfs = stats.binom.pmf(possible_values, n=n, p=p)

# Find the value(s) with the highest PMF
max_pmf = np.max(pmfs)
modes = possible_values[pmfs == max_pmf]

print(f"Binomial distribution mode(s): {modes}")

内容的提问来源于stack exchange，提问作者Adnan Tamimi