Numerical MLE in Python

1. From Maximization to Minimization

Let's recall the core idea of MLE. Given a set of i.i.d. data points $y_1,y_2,\dots ,y_n$ and a parametric model with probability density function (PDF) $f(y|\theta )$, the log-likelihood function is:

The Maximum Likelihood Estimate, ${\hat{\theta }}_{MLE}$, is the value of $\theta$ that maximizes this function:

Numerical optimization libraries are typically built to find the minimum of a function, not the maximum. Fortunately, any maximization problem can be trivially converted into a minimization problem. Maximizing a function $g(\theta )$ is equivalent to minimizing its negative, $-g(\theta )$.

Therefore, we can find the MLE by solving the following minimization problem:

This function, the negative log-likelihood, will be our "criterion function" for the optimizer.

2. A Practical Tool: scipy.optimize.minimize

The standard tool for numerical optimization in Python is the scipy.optimize.minimize function. It's a versatile function that provides access to several different optimization algorithms.

To get a practical feel for how it works, please watch the following video. It provides a clear, code-along tutorial for using minimize on both single-variable and multi-variable problems.

As you saw in the video, using scipy.optimize.minimize involves three main steps:

1. Define a criterion function: This Python function takes a vector of parameters as its first argument and returns a single scalar value to be minimized.
2. Provide an initial guess: The optimizer needs a starting point for its search for the optimal parameters.
3. Call minimize: Pass the function, the initial guess, and optionally other arguments like the data, bounds for the parameters, and constraints.

3. Implementing a Numerical MLE

Now, let's combine the theory of MLE with the practical tool of scipy.optimize.minimize. The process of implementing a numerical MLE can be broken down into clear steps.

Let's formalize the steps for our task:

1. Define the Model's PDF/PMF: Identify or define the probability distribution for your data.
2. Write the Log-Likelihood Function: Based on the PDF, write a Python function that computes the log-likelihood for a given set of parameters and data.
3. Create the Criterion Function: Create a wrapper function that takes the parameters and data, calls the log-likelihood function, and returns its negative value.
4. Optimize: Call scipy.optimize.minimize with your criterion function, an initial guess for the parameters, and your data passed via the args argument. It's also good practice to use the bounds argument to constrain parameters to their valid ranges (e.g., variance must be positive).
5. Extract the Results: The optimal parameters that constitute the MLE are in the .x attribute of the result object returned by minimize.

Hands-on: MLE for a Beta Distribution

Now it's time to put this into practice. We'll work with the Beta distribution, a continuous probability distribution defined on the interval $[0,1]$. It's governed by two positive shape parameters, $\alpha$ and $\beta$. Its PDF is:

where $B(\alpha ,\beta )$ is the Beta function. The log-likelihood for $n$ i.i.d. samples is:

Computing $\log B(\alpha ,\beta )$ can be done using log-gamma functions: $\log B(\alpha ,\beta )=\log \Gamma (\alpha )+\log \Gamma (\beta )-\log \Gamma (\alpha +\beta )$.

Your task is to write a Python script to find the MLE for $\alpha$ and $\beta$ for a set of data drawn from a Beta distribution.

Programming Exercise:

1. Import libraries: numpy, scipy.optimize, and scipy.special.
2. Generate data: Create a dataset of 1000 samples from a Beta distribution with known parameters, say $\alpha =2$ and $\beta =5$. Use numpy.random.beta.
3. Define the negative log-likelihood function:
   • This function should accept a list/array params (where params[0] is $\alpha$ and params[1] is $\beta$) and the data array.
   • Use scipy.special.loggamma to compute the log of the Beta function part.
   • Return the negative of the total log-likelihood.
4. Set initial guess and bounds:
   • Use an initial guess, e.g., [1, 1].
   • Since $\alpha >0$ and $\beta >0$, set bounds to be [(1e-6, None), (1e-6, None)].
5. Run the optimizer: Call scipy.optimize.minimize, passing the function, initial guess, data (via args), and bounds.
6. Print the results: Check the success flag and print the estimated parameters from result.x. Compare them to the true values of $\alpha =2$ and $\beta =5$.

import numpy as np
from scipy.optimize import minimize
from scipy.special import loggamma




# 1. True parameters
alpha_true = 2
beta_true = 5




# 2. Generate data
np.random.seed(42) # for reproducibility
data = np.random.beta(alpha_true, beta_true, size=1000)




# 3. Define the negative log-likelihood function
def neg_log_likelihood_beta(params, data):
    alpha, beta = params[0], params[1]
    



    # Log-likelihood components
    sum_log_x = np.sum(np.log(data))
    sum_log_1_minus_x = np.sum(np.log(1 - data))
    n = len(data)
    



    # Log of the Beta function
    log_beta_func = loggamma(alpha) + loggamma(beta) - loggamma(alpha + beta)
    



    # Total log-likelihood
    log_likelihood = (alpha - 1) * sum_log_x + \
                     (beta - 1) * sum_log_1_minus_x - \
                     n * log_beta_func
                     



    # Return the negative log-likelihood
    return -log_likelihood




# 4. Set initial guess and bounds
initial_guess = [1, 1]
bounds = [(1e-6, None), (1e-6, None)] # alpha > 0, beta > 0




# 5. Run the optimizer
result = minimize(
    neg_log_likelihood_beta,
    initial_guess,
    args=(data,),
    bounds=bounds,
    method='L-BFGS-B' # A good method for bound-constrained problems
)




# 6. Print the results
if result.success:
    estimated_params = result.x
    print(f"Optimization successful.")
    print(f"True parameters: alpha = {alpha_true}, beta = {beta_true}")
    print(f"Estimated MLE parameters: alpha = {estimated_params[0]:.4f}, beta = {estimated_params[1]:.4f}")
else:
    print("Optimization failed.")
    print(result.message)

Running this code, you should see that the estimated parameters are very close to the true values, demonstrating that our numerical MLE implementation works!

For Further Exploration: Gradient-Based Optimization

The minimize function is powerful, but for very complex models, its performance can be significantly improved by providing the optimizer with the analytical derivatives of the criterion function: the gradient (first derivative) and the Hessian (second derivative).

Algorithms like Newton-Raphson explicitly use these derivatives in their update rule. While calculating these by hand can be tedious and error-prone, modern automatic differentiation libraries like JAX can compute them for you.

Conclusion

In this lesson, you have acquired a powerful and versatile skill: finding Maximum Likelihood Estimates using numerical optimization. This frees you from the limitations of models with analytical solutions and allows you to fit virtually any custom statistical model to your data.

Key Takeaways:

• When an analytical solution for MLE is not feasible, we can use numerical optimization.
• The MLE problem is framed as a minimization problem by minimizing the negative log-likelihood.
• The scipy.optimize.minimize function is the primary tool for this task in Python.
• The implementation workflow involves defining a criterion function (the negative log-likelihood), providing an initial guess and bounds, and calling the optimizer.

Preview of the Next Lesson:

Now that we can find MLEs both analytically and numerically, we will turn our attention to one of their most elegant and useful features: the invariance property. This property provides a simple way to find the MLE for a function of a parameter if you already know the MLE of the parameter itself. For instance, if you have ${\hat{\sigma }}_{MLE}^2$, what is ${\hat{\sigma }}_{MLE}$? The invariance property gives a direct and intuitive answer, further solidifying our understanding of the MLE framework.