Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Maximum Likelihood Estimation › Implementation in NumPy and SciPy

Michael BrenndoerferOctober 25, 20267 min read

Machine Learning from ScratchAll chapters

Part 3 of 15Model Fundamentals and Diagnostics

Chapter 29Maximum Likelihood Estimation › Implementation in NumPy and SciPy

Reading section 7 of 10 · Implementation in NumPy and SciPy · 7 minSection 7 of 10 · 7 min

Section 7 of 10 · 7 min read

Maximum Likelihood Estimation › Implementation in NumPy and SciPySection view · Sections in this chapter7 / 10 ↓

Sections in this chapter10

Reading tools Definitions, links & notes

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Notes & comments

Notes & commentsLeave notes and join the discussionSelection toolbarShow comment actions when selecting text. Turn off for the browser menu only.Comment indicatorsShow small counts beside annotated passagesPassage highlightsGently shade annotated text selections

You will now implement the same calculations from primitives. First you compute the closed-form estimates directly. Then you compute a likelihood grid and run numerical optimization with SciPy. Finally, we compute curvature-based standard errors. The data are small and deterministic so that every printed value can be checked by hand.

The next cell imports SciPy's optimizer. The earlier hidden setup already imported NumPy, Matplotlib and the normal distribution and seeded the shared plot environment.

PythonIn [13]
CollapseExpandfrom scipy.optimize import minimize

You start with the closed-form estimates for the Bernoulli probability and the Gaussian mean and variance. These are the same formulas derived in the previous sections, now written as direct computations.

PythonIn [14]
CollapseExpandclicks = 12
visitors = 20
p_hat_closed = clicks / visitors

mu_hat_closed = y.mean()
sigma2_hat_closed = np.mean((y - mu_hat_closed) ** 2)
sigma_hat_closed = np.sqrt(sigma2_hat_closed)
sigma2_unbiased = np.var(y, ddof=1)

ConsoleOut [15]
CollapseExpandBernoulli closed-form MLE: p_hat = 0.6000
Gaussian closed-form MLE: mu_hat = 10.0000
Gaussian closed-form MLE: sigma2_hat = 0.0440 (sigma_hat = 0.2098)
Unbiased sample variance for comparison: 0.0550

Next we evaluate the Bernoulli log-likelihood on a grid and run a bounded numerical optimizer. For the Gaussian model, we optimize $\mu$ and $\log\sigma$ so that the standard deviation remains positive. Using $\log\sigma$ is a common practical trick. It turns a constrained optimization problem into an unconstrained one.

In the Gaussian optimizer below, the starting value for $\mu$ is set to the analytic sample mean mu_hat_closed.

PythonIn [16]
CollapseExpandp_grid = np.linspace(0.001, 0.999, 999)
bernoulli_ll_grid = clicks * np.log(p_grid) + (visitors - clicks) * np.log(
    1 - p_grid
)
p_hat_grid = p_grid[np.argmax(bernoulli_ll_grid)]


def bernoulli_neg_loglik(theta):
    p = theta[0]
    return -(clicks * np.log(p) + (visitors - clicks) * np.log(1 - p))


bernoulli_fit = minimize(
    bernoulli_neg_loglik,
    x0=np.array([0.5]),
    method="L-BFGS-B",
    bounds=[(1e-9, 1 - 1e-9)],
)
p_hat_numeric = bernoulli_fit.x[0]


def gaussian_neg_loglik(theta):
    mu, log_sigma = theta
    sigma = np.exp(log_sigma)
    return len(y) * (0.5 * np.log(2 * np.pi) + log_sigma) + np.sum(
        (y - mu) ** 2
    ) / (2 * sigma**2)


gaussian_fit = minimize(
    gaussian_neg_loglik,
    x0=np.array([mu_hat_closed, 0.0]),
    method="L-BFGS-B",
)
mu_hat_numeric, log_sigma_hat = gaussian_fit.x
sigma_hat_numeric = np.exp(log_sigma_hat)
sigma2_hat_numeric = sigma_hat_numeric**2

# Independent distribution fit and density implementation.
mu_reference, sigma_reference = norm.fit(y)
assert bernoulli_fit.success and gaussian_fit.success
assert np.isclose(p_hat_numeric, p_hat_closed, atol=1e-6, rtol=0)
assert np.allclose(
    [mu_hat_numeric, sigma_hat_numeric],
    [mu_reference, sigma_reference],
    atol=1e-5,
    rtol=0,
)
assert np.isclose(
    gaussian_neg_loglik(gaussian_fit.x),
    -norm.logpdf(y, loc=mu_hat_numeric, scale=sigma_hat_numeric).sum(),
)

ConsoleOut [17]
CollapseExpandBernoulli grid best p: 0.6000
Bernoulli numerical MLE: 0.6000
Gaussian numerical MLE: mu_hat = 10.0000, sigma2_hat = 0.0440
Optimizer success flags: Bernoulli=True, Gaussian=True
SciPy normal fit: mean=10.0000, variance=0.0440

The grid has spacing 0.001 and includes 0.6, so its answer matches p_hat_closed to four decimal places. The numerical optimizer should match the closed-form answer to several decimal places. If the optimizer reports failure, the estimate is not necessarily wrong, but it is a signal to inspect the starting value, scaling, and bounds.

Both routes recover the sample proportion and, for the Gaussian model, the optimizer recovers the sample mean and the maximum likelihood variance reported in the closed-form output above. Agreement among the grid search, the closed-form calculation, and the optimizer suggests the estimate is not an artifact of one particular numerical approach. When the routines disagree, the discrepancy usually points to grid resolution, an optimizer that stopped early, or a likelihood that is flat near its peak.

Finally, you compute approximate standard errors from the curvature formulas.

PythonIn [18]
CollapseExpand# Curvature-based standard errors from observed information
p_information = visitors / (p_hat_closed * (1 - p_hat_closed))
mu_information = len(y) / sigma2_hat_closed
sigma2_information = len(y) / (2 * sigma2_hat_closed**2)

Each information value summarizes how sharply the log-likelihood peaks at the estimate. Larger information implies a narrower peak and a smaller approximate standard error. The standard errors below are the square roots of the reciprocals of these information values. For the Bernoulli probability, the information is $n/[\hat p(1-\hat p)]$ . For the Gaussian mean, the information is $n/\hat\sigma^2$ . For the Gaussian variance, the asymptotic standard error is $\hat\sigma^2\sqrt{2/n}$ , which follows from the observed information $n/(2\hat\sigma^4)$ and the square root of its reciprocal. The two are related by $\mathrm{se}(\hat\sigma^2) \approx \hat\sigma^2\sqrt{2/n} = \sqrt{2\hat\sigma^4/n}$ , and that approximate value is computed from the observed information in the code below rather than hardcoded. Explicitly,

\sqrt{\frac{2\hat\sigma^4}{n}} = \sqrt{\hat\sigma^4}\sqrt{\frac{2}{n}} = \hat\sigma^2\sqrt{\frac{2}{n}}.

where:

$\hat\sigma^2$ : the maximum likelihood variance estimate
$n$ : the number of observations
$\frac{n}{2\hat\sigma^4}$ : the observed information for the Gaussian variance

PythonIn [19]
CollapseExpandp_se_closed = np.sqrt(1 / p_information)
mu_se_closed = np.sqrt(1 / mu_information)
sigma2_se_closed = np.sqrt(1 / sigma2_information)

ConsoleOut [20]
CollapseExpandApproximate standard error for p_hat: 0.1095
Approximate standard error for mu_hat: 0.0938
Approximate standard error for sigma2_hat: 0.0278
Root of observed information for p_hat: 9.1287

The Bernoulli standard error is the same value we derived by hand from clicks and visitors, and it reflects the moderate amount of information in twenty visits. The Gaussian mean standard error is small in absolute terms, but with only five observations it relies on the asymptotic approximation. Here the variance standard error is large relative to the variance estimate. Comparing this with the mean estimate is not a general ranking of how difficult means and variances are to estimate: their scales and units differ. Variance intervals also need care in small samples.

Key ParametersLink Copied

The key parameters that govern the implementation are:

clicks and visitors: the number of successes and trials for the Bernoulli example. The closed-form estimate is clicks / visitors; the ratio is only meaningful when visitors is positive.
p_grid: the sequence of candidate probabilities between $0.001$ and $0.999$ used by the grid search. A finer grid can reduce discretization error but increases cost, and the bounds keep $p$ away from the exact boundary where the logarithm is undefined.
Optimizer starting value for p: 0.5 for the Bernoulli fit; a plausible interior value helps the quasi-Newton method converge.
y: the observed Gaussian vector. The closed-form mean is y.mean(), and the maximum likelihood variance is the average squared deviation from that mean, np.mean((y - mu_hat_closed)**2).
Optimizer starting value x0 = [mu_hat_closed, 0.0]: the first entry is an initial guess for $\mu$ (the sample mean) and the second is an initial guess for $\log\sigma$ . The log parameterization keeps the standard deviation positive.
method="L-BFGS-B": a bounded quasi-Newton optimizer that works well for smooth lower-dimensional problems.
Observed information: visitors / (p_hat * (1 - p_hat)) for the Bernoulli probability, len(y) / sigma2_hat for the Gaussian mean, and len(y) / (2 * sigma2_hat**2) for the Gaussian variance. The information formulas are exact for the stated interior fits; using inverse information to estimate sampling uncertainty is the large-sample approximation.

Key MethodsLink Copied

np.mean computes the sample mean by summing the observations and dividing by $n$ . np.var(y, ddof=0) would compute the maximum likelihood variance, while np.var(y, ddof=1) computes the unbiased sample variance. The difference is the denominator $n$ versus $n-1$ . np.argmax returns the index of the largest value in an array, which is how you select the best grid point.

scipy.optimize.minimize minimizes a scalar function of one or more variables. We pass the negative log-likelihood because the optimizer minimizes. The method="L-BFGS-B" option supports bounds. In these calls, no jac is supplied, so SciPy estimates gradients numerically; callers can also supply an analytic gradient. norm.logpdf evaluates normal log-densities, whose sum supplies a reference for our primitive objective. norm.fit(y) estimates the normal location and scale, which here are the sample mean and maximum likelihood standard deviation.

The optimizer returns an object with x, the estimated parameters, and success, a Boolean flag indicating whether the algorithm believes it converged. Convergence flags are useful but not sufficient. For an unconstrained interior fit, check that the gradient is near zero and inspect the Hessian of the negative log-likelihood. Positive definiteness is a sufficient local minimum check, not proof of global optimality; bounded solutions instead require the corresponding constrained conditions. Compare reasonable starting values as well.

Comments

Add to the conversation

A question, a correction, or an example from your own work.

No comments yet. Be the first to share your thoughts!

Key ParametersLink Copied

Key MethodsLink Copied

Comments

Add to the conversation

Stay up to date