Similar to How to get the mode of distribution in scipy.stats, The scipy.stats library has functions to find the mean and median of a fitted distribution but not mode. If I have the parameters of a distribution after fitting to data, how can I find the mode of the fitted distribution?

I would like to know how I would get the mode of a log-normal distribution in Scipy after the update? I've not been able to understand.

This is in reference to:

Currently, this is the only solution I have

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# Define the parameters for the underlying normal distribution
mu, sigma = 0, 1  # Mean and standard deviation of the normal distribution
# Theoretical mode based on the formula for a log-normal distribution
theoretical_mode = np.exp(mu - sigma**2)

lognorm_dist = stats.lognorm(sigma, scale=np.exp(mu))

x = np.linspace(0, 10, 1_000) # Create an array of x values to evaluate the PDF
y = lognorm_dist.pdf(x) # Get the PDF values for each x

# Find the mode (maximum of the PDF)
mode = x[np.argmax(y)]  # The x value corresponding to the maximum PDF

print(f"Numerical Mode of the log-normal distribution: {mode}")
print(f"Theoretical Mode of the log-normal distribution: {theoretical_mode}")

# Plot
plt.plot(x, y, label='Log-Normal Distribution')
plt.axvline(mode, color='red', linestyle='--', label=f'Numerical Mode: {mode}')
plt.axvline(theoretical_mode, color='green', linestyle='--', label=f'Theoretical Mode: {theoretical_mode}')
plt.legend()
plt.show()

Similar to How to get the mode of distribution in scipy.stats, The scipy.stats library has functions to find the mean and median of a fitted distribution but not mode. If I have the parameters of a distribution after fitting to data, how can I find the mode of the fitted distribution?

I would like to know how I would get the mode of a log-normal distribution in Scipy after the https://github.com/scipy/scipy/pull/21050 update? I've not been able to understand.

This is in reference to: https://github.com/scipy/scipy/issues/14895#issuecomment-1186703157

Currently, this is the only solution I have

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# Define the parameters for the underlying normal distribution
mu, sigma = 0, 1  # Mean and standard deviation of the normal distribution
# Theoretical mode based on the formula for a log-normal distribution
theoretical_mode = np.exp(mu - sigma**2)

lognorm_dist = stats.lognorm(sigma, scale=np.exp(mu))

x = np.linspace(0, 10, 1_000) # Create an array of x values to evaluate the PDF
y = lognorm_dist.pdf(x) # Get the PDF values for each x

# Find the mode (maximum of the PDF)
mode = x[np.argmax(y)]  # The x value corresponding to the maximum PDF

print(f"Numerical Mode of the log-normal distribution: {mode}")
print(f"Theoretical Mode of the log-normal distribution: {theoretical_mode}")

# Plot
plt.plot(x, y, label='Log-Normal Distribution')
plt.axvline(mode, color='red', linestyle='--', label=f'Numerical Mode: {mode}')
plt.axvline(theoretical_mode, color='green', linestyle='--', label=f'Theoretical Mode: {theoretical_mode}')
plt.legend()
plt.show()

• scipy

1 Answer 1

As of SciPy 1.15, there is a new random variable infrastructure, and random variable objects have a mode method.

One way to get the mode of a lognormal distribution is to instantate an underlying normal random variable, exponentiate it (since the log of the resulting random variable is normal), and find the mode of the resulting random variable.

from scipy import stats
import matplotlib.pyplot as plt

mu, sigma = 0, 1
Z = stats.Normal(mu=mu, sigma=sigma)
Y = stats.exp(Z)
mode = Y.mode()  # 0.3678794390188982

Y.plot(t=('x', 0, 5))
plt.plot(mode, Y.pdf(mode), '*')

If you already have your lognormal distribution parameterized in terms of the s shape parameter and scale parameter of the old stats.lognorm distribution, you might prefer to expose stats.lognorm with the new infrastructure using make_distribution.

import numpy as np
from scipy import stats

# expose `stats.lognorm` with new infrastructure
LogNormal = stats.make_distribution(stats.lognorm)

# lognorm_dist = stats.lognorm(sigma, scale=np.exp(mu))
X = LogNormal(s=sigma) * np.exp(mu)
X.mode()  # 0.36787943662147315

Both of these use numerical optimization to perform the calculation, so they are not as precise as your theoretical value. However, that's an inherent limitation of direct numerical optimization rather than a shortcoming of the implementation; the PDFs at these points are indistinguishable from one another in double precision.

theoretical_mode = np.exp(mu - sigma**2)
(Y.pdf(Y.mode()) 
 == X.pdf(X.mode()) 
 == stats.lognorm(sigma, scale=np.exp(mu)).pdf(theoretical_mode))
# True, at least on my machine

For more information, see the Random Variable Transition Guide and method documentation of mode.