• P(y | θ): Probability of data y when parameter θ is known. θ fixed, y unknown.
• L(θ | y): Likelihood of parameter θ when data y are observed. y fixed, θ unknown.

Same mathematical formula — completely different interpretation.

Likelihood does not sum to 1 and is not a probability over θ. For continuous distributions L is the density value — which can exceed 1 (set σ very small in Stage 2 to see the density curve rise above 1).

The likelihood of n observations is the product of individual densities: L(θ|y₁…yₙ) = ∏ᵢ f(yᵢ|θ)

With many small values (e.g. 0.04 × 0.09 × … × 0.06) the product becomes extremely small — numerical underflow. Stage 2 shows this live: the product of the red lines is displayed in the info text and shrinks rapidly.

The logarithm converts the product into a sum: ℓ(θ) = Σᵢ log f(yᵢ|θ) — it is a monotonically increasing function, so the maximum lies at exactly the same point θ whether you maximize L(θ) or ℓ(θ).

The key practical advantage: instead of multiplying thousands of tiny numbers together (numerical underflow — the computer simply returns 0), you now add manageable negative numbers, with no loss of precision. MLE therefore always maximizes ℓ(θ) = log L(θ).

Maximum Likelihood Estimation finds θ that makes the observed data most plausible: θ = argmax ℓ(θ|y)

This is not a verdict about θ itself — only about its compatibility with the data. Other θ values are not impossible, just less plausible.

Analytical MLE estimators (Normal): μ = y    σ² = Σ(yᵢ−y)²/n

Why is σ² biased? The estimator centers on y itself, and y is, by construction, exactly the value that minimizes Σ(yᵢ−c)² over all possible c (the least-squares property of the mean). That makes Σ(yᵢ−y)² smaller on average than Σ(yᵢ−μ)² with the true μ. Exactly: E[σ²] = (n−1)/n · σ² — the effect is strongest for small n. The unbiased estimator s² = Σ(yᵢ−y)²/(n−1) corrects for this (Bessel's correction) — that's why R's var() divides by n−1.

In Stage 1: μ = y (one point). In Stage 2: μ = y (all points).

MLE can estimate multiple parameters simultaneously — a key advantage over simple moment estimators.

For the Normal distribution there are two unknown parameters: μ (location) and σ (spread). The log-likelihood becomes a surface over the (μ, σ)-space — a mountain landscape instead of a curve. The joint MLE lies at the peak: μ = y    σ = √(Σ(yᵢ−y)²/n)

In Stage 2: activate "Vary σ too" → the heatmap shows the 2D landscape, the bright spot is the joint peak. This principle scales to arbitrarily many parameters — this is how lm() and glm() in R work internally.

MLE works for any parametric family — not just Normal. This is the foundation of GLMs:

• Poisson (training hours/week): λ = y  |  Log-likelihood: Σᵢ [yᵢ·log(λ) − λ]
• Bernoulli (min. performance 0/1): p = proportion of ones  |  Log-likelihood: Σᵢ [yᵢ·log(p) + (1−yᵢ)·log(1−p)]

Switch between families in Stage 3: the likelihood landscape changes its shape — parabola for Bernoulli, asymmetric for Poisson — but the mechanism is identical: find the peak.

The maximized log-likelihood value ℓ̂ can be used directly for model comparison:

• AIC = −2·ℓ̂ + 2k  (k = number of parameters)
• BIC = −2·ℓ̂ + k·log(n)

Smaller AIC/BIC = better fit at equal complexity. AIC penalizes less strongly than BIC — with large n, BIC favors more parsimonious models.

Important caveat: AIC/BIC are directly comparable only between models of the same distributional class — e.g. several count-data models (Poisson vs. Negative-Binomial) or several normal models. Between fundamentally different families (e.g. Bernoulli vs. Normal), the values are not comparable without further adjustment, since they rest on different data scales/likelihood normalizations.

Application example: Does Poisson or Negative-Binomial fit the training-hours data better? Same data, same distributional class (count data) — here AIC/BIC decide. In Stage 3, both are computed live.