Estimation, likelihood, and Bayesian updating

Fitting a parameter requires more than finding a large number on a curve. We must say which data mechanism the curve represents, which quantities are unknown, and what loss or uncertainty interpretation accompanies the reported result. This module derives Bernoulli and Gaussian maximum likelihood, separates bias from variance, performs a fully normalised beta–Bernoulli update, and traces regularisation coefficients back to noise and prior assumptions. A likelihood, a sampling distribution and a posterior are three different mathematical objects.

Retrieval check: use Module 16 for derivatives, Module 17 for densities and integrals, Module 23 for moments, and Module 24 for sampling limits. All labs use only the Python standard library. The continuous parameter and Gaussian examples make calculus prerequisites explicit even when the observed data are discrete.

A statistical model is a collection of possible probability laws P_θ indexed by a parameter θ in a stated parameter space Θ. The data are generated under one of these laws if the model is correctly specified. For Bernoulli observations, Θ=[0,1] and θ=p is a success rate. For a Gaussian location model with known positive variance, Θ is the real line and θ=μ. If the variance is also unknown, the parameter is a pair (μ,σ²) with σ²>0, and factors that were constant in μ may now depend on the parameter.

In a frequentist sampling analysis θ is fixed but unknown, while the sample and any statistic computed from it are random. An estimator T(X₁,…,X_n) is a specified rule before data are observed. Its realised numerical output t is an estimate. The estimator has a sampling distribution indexed by θ; the estimate is one outcome. Confusing these levels turns a calculation such as E_θ[T]=θ into a claim that every observed estimate equals the truth, which unbiasedness never says.

In Bayesian analysis a prior distribution is assigned to the unknown parameter and updated by observed data through a joint model. Conditional on the parameter, the observations may be independent; after integrating over a common uncertain parameter, they can become dependent. State which conditioning level an independence assertion concerns. A Bayesian posterior mean is still a data-dependent estimator and can be assessed by frequentist sampling risk at each fixed θ, without equating that risk with a posterior probability statement.

Identifiability means different parameter values imply different observable laws. If P_θ=P_θ′ while θ≠θ′, no number of observations from that law can identify which label was used without extra restrictions. Consider X distributed as Gaussian with known variance and mean θ² for θ in the real line. Values θ and −θ give identical laws, so the signed θ is unidentifiable. Restricting Θ to nonnegative values makes the mean-to-parameter map injective, though uncertainty near zero still requires care.

The sampling scheme is part of the model. IID Bernoulli likelihood assumes a common rate and conditional independence of the recorded trials. Duplicated rows or changing rates do not obey that product likelihood merely because each row contains zero or one. Selection depending on unobserved outcomes can also alter the likelihood of the observed sample. If a model is an approximation to reality, its fitted parameter may describe the best fit within that approximation rather than a literal population mechanism. Later regression lessons assess residual and calibration evidence for such mismatches.

A useful estimation report specifies the observable variables, sample unit, parameter space, dependence assumptions, estimator rule and target. These statements decide whether a boundary is admissible, whether a derivative equation suffices, and which expectation defines bias. A parameter value outside its space is invalid even if an unconstrained optimisation formula returns it. Conversely a legal boundary optimum must not be discarded just because no interior derivative is zero there.

For observed data x, likelihood is L(θ;x)=p_θ(x), using a mass function for discrete observations or a density with respect to a common reference measure for continuous observations. The numerical expression may be the same as a probability-model expression, but likelihood fixes x and varies θ. A probability law instead fixes θ and sums or integrates over possible x. A likelihood does not automatically integrate to one over Θ, and normalising it over a parameter coordinate is itself an extra modelling choice.

For conditionally independent observations, L(θ;x)=∏p_θ(x_i). The log-likelihood ℓ=Σlog p_θ(x_i) has the same maximisers whenever the likelihood is positive because log is strictly increasing. Impossible observations at a parameter value give L=0 and ℓ=−∞. Logarithms transform products to sums and prevent many underflow failures, but do not repair a wrong independence assumption. Use the extended-value boundary convention explicitly rather than evaluate expressions such as zero times log zero as an accidental floating-point NaN.

For an ordered Bernoulli sample with k successes and n−k failures, L(p)=p^k(1−p)^(n−k). Observing only the count K=k gives the binomial likelihood binom(n,k)p^k(1−p)^(n−k). The combinatorial factor is constant in p for fixed n,k, so both have the same p maximiser and posterior shape under the same prior. Their data probabilities differ because the two observations describe different events. A factor can cancel in one parameter-fitting task and still matter when computing probabilities or comparing other models.

Dropping constants is legitimate only after naming the argument being optimised or normalised. In a Gaussian mean fit with known σ, the prefactor is constant in μ. If σ is unknown, the term σ^(−n) is parameter-dependent and cannot be discarded. In a posterior update, a factor independent of θ cancels between numerator and evidence, but the posterior must still have a valid normalising integral. In model comparison, factors depending on the competing model can change the result even if each was irrelevant to its own within-model optimiser.

For a finite candidate list, relative weights can be computed with log-sum-exp: log Σ exp a_j = m+log Σ exp(a_j−m), where m=max a_j is finite. This stabilises arithmetic. Turning those weights into a posterior over candidates still requires prior masses and a specified candidate model; the numerical normalisation alone does not supply a prior interpretation. Continuous evidence is an integral, so a sum over an arbitrary grid additionally needs quadrature weights and a justified approximation.

Likelihood values for continuous data are densities, not probabilities of observing an exact real-valued point. Densities can exceed one and change numerical units under transformations of observations. An invertible data transformation contributes its appropriate Jacobian; if that Jacobian is independent of θ, relative fitting is preserved. Transforming a parameter is a different operation: a likelihood is composed with the coordinate map, whereas a prior or posterior density transforms with a Jacobian. Section 6 uses this distinction to explain why continuous MAP modes depend on coordinates.

For Bernoulli data with 0<k<n, ℓ(p)=k log p+(n−k)log(1−p) on 0<p<1. Its derivative is k/p−(n−k)/(1−p). Setting this to zero gives k(1−p)=(n−k)p and p_hat=k/n. The second derivative is −k/p²−(n−k)/(1−p)²<0, establishing strict concavity and the unique interior maximum. Both boundaries give zero likelihood when both outcomes occur. A complete maximisation argument therefore checks curvature and admissible boundaries, rather than treating any derivative root as sufficient.

When k=0 and n>0, L(p)=(1−p)^n is maximised at p=0 on [0,1]; when k=n it is maximised at p=1. On an open parameter space (0,1), these data instead have a supremum approached at a missing endpoint and no maximiser. With no data, n=0, every parameter has likelihood one and the data supply no unique estimate. The parameter-space choice changes whether a maximum exists, even though the algebraic likelihood expression is unchanged.

For IID Gaussian observations with known variance σ²>0, the full log-likelihood is −n log(√(2π)σ)−Σ(x_i−μ)²/(2σ²). Maximising over μ minimises the residual sum of squares. Expand around x_bar: Σ(x_i−μ)²=Σ(x_i−x_bar)²+n(μ−x_bar)², because the centred residuals sum to zero. The second term is nonnegative and vanishes only at μ=x_bar for n>0. This identity proves the mean MLE globally, with derivative n(x_bar−μ)/σ² giving the same result.

If variance v=σ² is also unknown and the centred residual sum R is positive, substitute μ=x_bar and maximise −(n/2)log v−R/(2v) plus constants. The derivative is −n/(2v)+R/(2v²), hence v_hat=R/n. The log-likelihood increases below R/n and decreases above it, and tends to −∞ at both variance extremes, so this is the global maximum. The MLE denominator is n; an unbiased variance estimator uses n−1 for a different purpose, derived next.

If all Gaussian observations are identical, R=0. Setting μ equal to that value and letting v tend to zero makes the likelihood unbounded above on v>0. There is no finite positive-variance MLE, and substituting v_hat=0 would leave the stated parameter space. A numerical variance floor changes the optimisation problem; it can be a practical model restriction but must be disclosed. Likewise a unique MLE in these elementary models does not justify claiming uniqueness or existence for every mixture, nonlinear or high-dimensional model.

Least squares emerges here from a Gaussian observation model. Other noise models give other objectives. IID Laplace location noise with known scale b>0 has likelihood proportional to exp(−Σ|x_i−μ|/b), so a sample median minimises absolute deviations. With an even number of observations an entire interval of medians may minimise that loss. Model choice, target and optimisation geometry therefore interact; selecting a loss is a statistical assumption as well as a computational choice.

For an estimator T of a fixed scalar θ with finite second moment, bias is E_θ[T]−θ, variance is E_θ[(T−E_θT)²], and mean squared error is E_θ[(T−θ)²]. Write T−θ=(T−E_θT)+(E_θT−θ), square and take expectations. The cross term is zero, giving MSE_θ(T)=Var_θ(T)+Bias_θ(T)². These are repeated-sample properties at a stated parameter, not quantities determined by one realised absolute error. The subscript reminds us that a comparison can change with the true θ.

For IID Bernoulli data, T=K/n is unbiased and has variance p(1−p)/n. Its MSE equals that variance. For a fixed Beta(a,b) posterior mean T_B=(K+a)/(n+a+b), frequentist expectation is (np+a)/(n+a+b), so bias is [a−(a+b)p]/(n+a+b), while variance is np(1−p)/(n+a+b)². Shrinkage reduces this variance but introduces parameter-dependent bias. Whether total MSE improves follows from the full sum, not from the variance alone or the word “Bayesian.”

For IID finite-variance X_i, let R=Σ(X_i−x_bar)². The identity R=Σ(X_i−μ)²−n(x_bar−μ)² follows by expanding squares. Expectations give E[R]=nσ²−n(σ²/n)=(n−1)σ². Thus S²=R/(n−1) is unbiased for n>1. The derivation uses a common mean, variance and zero independent covariances; treating correlated rows as independent can invalidate the denominator correction. Gaussianity is unnecessary for this expectation identity, although particular exact sampling distributions require it.

For x=(1,2,4), x_bar=7/3 and R=14/3, giving Gaussian variance MLE 14/9 and unbiased sample variance 7/3. The latter is larger by n/(n−1). “Unbiased” refers to estimating variance, not to its square root: taking a nonlinear square root generally changes bias. An n−1 correction also does not guarantee smaller MSE than every biased alternative. The criterion being optimised must be named before calling an estimator better.

Consistency means T_n converges in probability to θ as sample size grows. The Bernoulli proportion is consistent by Module 24’s weak law. With fixed positive a,b, T_B−K/n=[a−(a+b)K/n]/(n+a+b) tends uniformly to zero because K/n stays in [0,1], so its posterior-mean rule is also consistent. Finite-sample bias and consistency can coexist. A fixed prior loses relative weight here; a prior strength increasing proportionally with n can preserve nonzero limiting shrinkage and need not give the same consistency result.

Unbiasedness alone also does not imply consistency: an estimator that always returns the first Bernoulli observation is unbiased for p but keeps variance p(1−p) as n grows. Conversely a small vanishing bias with vanishing variance gives consistency by a suitable squared-error bound. Do not transfer consistency or unbiasedness from a sample mean to arbitrary maximum-likelihood algorithms without a proof and model conditions. General MLE asymptotics involve identifiability, regularity and possible boundary complications beyond the elementary examples established here.

Risk comparisons should also hold the data budget and estimator rule fixed. Selecting a rule after inspecting the same outcomes creates a new data-dependent estimator whose expectation is not automatically the formula for either original rule. Lab 2 specifies both rules and both synthetic rates before generating datasets; model selection protocols in the next lesson make that distinction operational.

A Beta(a,b) prior for 0<p<1 has density π(p)=p^(a−1)(1−p)^(b−1)/B(a,b), with a,b>0 and B(a,b)=∫_0^1 p^(a−1)(1−p)^(b−1)dp. The positive parameters make the endpoint singularities integrable if present. For positive integer a,b, repeated integration by parts yields B(a,b)=(a−1)!(b−1)!/(a+b−1)!. The prior is a distribution over the rate, not extra observed trials; its shape expresses a modelling choice made before the claimed evidence update.

Multiply the Bernoulli likelihood by the prior density. Powers add, giving p^(a+k−1)(1−p)^(b+n−k−1), whose integral is B(a+k,b+n−k). Bayes therefore yields the properly normalised Beta(a+k,b+n−k) posterior. For the ordered sample, evidence is B(a+k,b+n−k)/B(a,b); for the count it also includes binom(n,k). That data-constant factor cancels from the posterior, but remains in the probability of the observed count. Evidence must be positive and finite for this conditioning step.

Using the ratio B(A+1,B)/B(A,B) gives posterior mean A/(A+B). The second-moment ratio gives A(A+1)/[(A+B)(A+B+1)], so posterior variance is AB/[(A+B)²(A+B+1)], where A=a+k and B=b+n−k. This is conditional uncertainty over p after these data under the specified prior and model. It is different from the frequentist sampling variance of the rule (K+a)/(n+a+b), which averages over new datasets at fixed p.

For a future Bernoulli trial conditionally independent of the recorded trials given p, posterior predictive success probability is ∫p π(p|x)dp=E[p|x]=A/(A+B). This prediction integrates parameter uncertainty rather than plugging in a density mode. Future trials are conditionally independent given p but generally dependent under the posterior predictive mixture: their shared uncertain rate gives covariance Var(p|x)>0. Thus the multiple-future-trial count is not automatically binomial at the posterior mean; a beta-binomial mixture captures that additional variation.

An equal-tail 95% credible interval uses posterior CDF quantiles .025 and .975. It contains posterior probability .95 under the chosen prior and likelihood; endpoints in Lab 1 come from numerical inversion checked by CDF mass. This does not automatically guarantee 95% frequentist coverage at every fixed p. Module 26 distinguishes repeated-sampling confidence intervals. A uniform prior is uniform in a specific coordinate, and credible conclusions can be sensitive to prior or likelihood misspecification, especially with limited data.

The predictive distinction can be quantified for r future trials. Let S be their success count and m=A/(A+B). Conditional on p, its mean is rp and variance rp(1−p). Total variance therefore gives Var(S|x)=rE[p(1−p)|x]+r²Var(p|x). Since E[p(1−p)|x]=m(1−m)−Var(p|x), this becomes r m(1−m)+r(r−1)Var(p|x). The first term is the plug-in binomial variance; the second is nonnegative additional uncertainty from sharing an unknown rate. It vanishes for one future trial, explaining why the single predictive probability is simple even though the joint future law differs. Under Beta(10,4), ten future trials have mean 50/7 and variance 10(5/7)(2/7)+90(2/147), about 3.265306 rather than the plug-in 2.040816. This is a conditional prediction calculation, not variation of the posterior mean across newly sampled training datasets.

The explorer uses integer prior parameters so its beta density and CDF can be computed transparently without an external statistical library. It plots relative likelihood in one panel and a normalised posterior density in the other, with separate labelled vertical scales. A relative-likelihood height of one is not a posterior probability. MLE, posterior mean, density mode and credible interval are labelled separately, including all-success, all-failure and uniform-posterior cases.

Sequential updating with genuinely new conditional-independent batches gives the same posterior as one combined update because their counts add. Updating with the same batch twice is different: it squares the likelihood and invents repeated evidence. For the eight/two data, reusing them after Beta(10,4) falsely produces Beta(18,6), not a justified refinement from new data. A prior estimated from the same data also requires a disclosure and analysis appropriate to that construction; it cannot silently be presented as fixed independent prior information.

Maximum a posteriori estimation chooses a mode of a posterior mass function or of a density in a stated continuous coordinate. For a density, maximising L(θ)π(θ) is equivalent to minimising negative log-likelihood plus negative log-prior, up to constants independent of θ. A density mode is not the point with positive posterior probability in a continuous model: individual points have probability zero. With multimodal or flat densities, a MAP summary can be nonunique and need not resemble a posterior mean or a useful predictive summary.

For Beta(A,B) with A,B>1, differentiate (A−1)log p+(B−1)log(1−p) to obtain mode (A−1)/(A+B−2). If A=1 and B>1, the maximum-density endpoint is zero in the continuous extension to [0,1]; if B=1 and A>1 it is one. For A=B=1 the density is uniform, so every point ties. With positive parameters below one, endpoint densities can diverge and require a careful mode convention. Our widget restricts integer shapes at least one to avoid hiding these singular cases.

Consider IID observations y_i∼N(μ,σ²) with known σ²>0 and prior μ∼N(0,τ²), τ²>0. Negative log posterior is Σ(y_i−μ)²/(2σ²)+μ²/(2τ²) plus a μ-constant. Multiplying by 2σ² gives Σ(y_i−μ)²+λμ² with λ=σ²/τ². Differentiation gives μ_MAP=n y_bar/(n+λ). If the data term is average squared error instead, the matching coefficient is λ/n. Changing sum to average while retaining the same penalty coefficient changes the prior-equivalent objective.

For a vector coefficient with independent zero-mean Gaussian prior components, the negative log prior is proportional to its squared Euclidean norm. Independent Laplace prior components instead produce an absolute-value penalty. This connection needs both the stated coefficient prior and observation model; a Gaussian noise likelihood supplies squared residuals, while Laplace noise supplies absolute residuals. Penalising an intercept, changing feature units or changing a noise scale also changes the probabilistic interpretation. Later regression and optimisation modules treat design matrices and constrained algorithms, but the scalar derivation already fixes the coefficient convention.

Continuous MAP depends on parameter coordinates. If p has Beta(A,B) posterior and η=log[p/(1−p)], the η-density includes the Jacobian dp/dη=p(1−p). Its shape as a function of p is therefore proportional to p^A(1−p)^B, with η-mode corresponding to p=A/(A+B), whereas the p-density mode is (A−1)/(A+B−2) when interior. The posterior distribution and properly transformed event probabilities remain consistent; density heights and their mode are coordinate-dependent. A prior should not be called neutral solely because it is flat in one convenient coordinate.

MAP regularisation is a modelling interpretation, not a guarantee of calibrated uncertainty or generalisation. A point optimiser discards much of the posterior, and a penalty tuned from validation outcomes must be described with that selection protocol. If σ is fitted too, its normalising term must return to the objective. If the prior depends on θ through an unknown hyperparameter, its normalisation can matter for fitting that hyperparameter. Keeping only familiar squared or absolute terms without checking which arguments vary can silently change the intended model.

Lab 1 · Likelihood and a normalised beta update

Download lab1_likelihood_and_beta_update.py

"""CPU-only: likelihood estimates and an integer-parameter beta posterior."""
import math


def beta_cdf(x, a, b):
    """For positive integer a,b, I_x(a,b) is a binomial upper tail."""
    if x <= 0:
        return 0.0
    if x >= 1:
        return 1.0
    n = a + b - 1
    return math.fsum(math.comb(n, j) * x**j * (1-x)**(n-j)
                     for j in range(a, n + 1))


def beta_quantile(u, a, b):
    lo, hi = 0.0, 1.0
    for _ in range(70):
        mid = (lo + hi) / 2
        if beta_cdf(mid, a, b) < u:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


if __name__ == "__main__":
    n, k, a, b = 10, 8, 2, 2
    post_a, post_b = a+k, b+n-k
    mle, mean = k/n, post_a/(post_a+post_b)
    mode = (post_a-1)/(post_a+post_b-2)
    print("Ordered Bernoulli likelihood: p^8(1-p)^2; count factor cancels for p fitting")
    for p in (.2, .5, .7, .8, .9):
        print(f"p={p:.1f}: likelihood={p**k*(1-p)**(n-k):.9f}")
    low, high = (beta_quantile(u, post_a, post_b) for u in (.025, .975))
    assert abs(beta_cdf(high, post_a, post_b)-beta_cdf(low, post_a, post_b)-.95) < 1e-12
    print(f"MLE={mle:.6f}; posterior Beta({post_a},{post_b})")
    print(f"posterior mean / next-trial predictive={mean:.6f}; posterior mode={mode:.6f}")
    print(f"equal-tail 95% credible interval=[{low:.6f}, {high:.6f}] (numerical inversion)")
    xs = [1., 2., 4.]
    mu = math.fsum(xs)/len(xs)
    rss = math.fsum((x-mu)**2 for x in xs)
    print(f"Gaussian data {xs}: mean MLE={mu:.6f}; variance MLE={rss/3:.6f}; unbiased variance={rss/2:.6f}")
    print("All failures: Bernoulli MLE=0; Beta(2,12) mean=1/7, mode=1/12")

Ordered Bernoulli likelihood: p^8(1-p)^2; count factor cancels for p fitting
p=0.2: likelihood=0.000001638
p=0.5: likelihood=0.000976562
p=0.7: likelihood=0.005188321
p=0.8: likelihood=0.006710886
p=0.9: likelihood=0.004304672
MLE=0.800000; posterior Beta(10,4)
posterior mean / next-trial predictive=0.714286; posterior mode=0.750000
equal-tail 95% credible interval=[0.461868, 0.909080] (numerical inversion)
Gaussian data [1.0, 2.0, 4.0]: mean MLE=2.333333; variance MLE=1.555556; unbiased variance=2.333333
All failures: Bernoulli MLE=0; Beta(2,12) mean=1/7, mode=1/12

Derive the three Bernoulli summaries before running. Check the interval’s posterior mass, then explain the two Gaussian variance denominators and the all-failure boundary.

Lab 2 · Repeated-sample bias variance and MSE

Download lab2_bias_variance_and_mse.py

"""Repeated independent datasets compare specified frequentist risks."""
import math
import random


def summarize(values, truth):
    avg = math.fsum(values)/len(values)
    var = math.fsum((v-avg)**2 for v in values)/len(values)
    mse = math.fsum((v-truth)**2 for v in values)/len(values)
    assert abs(mse-var-(avg-truth)**2) < 1e-12
    return avg-truth, var, mse


if __name__ == "__main__":
    rng = random.Random(25025)
    n, replications = 10, 6000
    print(f"seed=25025; independent synthetic datasets={replications}; n={n}; Beta(2,2) fixed before draws")
    for p in (.3, .95):
        counts = [sum(rng.random() < p for _ in range(n)) for _ in range(replications)]
        estimates = {"MLE": [k/n for k in counts], "posterior mean": [(k+2)/(n+4) for k in counts]}
        for name, values in estimates.items():
            bias, var, mse = summarize(values, p)
            model_bias = 0. if name == "MLE" else (2-4*p)/(n+4)
            model_var = p*(1-p)/n if name == "MLE" else n*p*(1-p)/(n+4)**2
            print(f"p={p:.2f}, {name}: empirical bias={bias:.6f}, variance={var:.6f}, MSE={mse:.6f}")
            print(f"  exact sampling bias={model_bias:.6f}, variance={model_var:.6f}, MSE={model_var+model_bias**2:.6f}")
    print("A posterior mean can be assessed by repeated-sample risk; its Bayesian interpretation does not guarantee lower risk at every fixed p.")

seed=25025; independent synthetic datasets=6000; n=10; Beta(2,2) fixed before draws
p=0.30, MLE: empirical bias=0.000417, variance=0.021185, MSE=0.021185
  exact sampling bias=0.000000, variance=0.021000, MSE=0.021000
p=0.30, posterior mean: empirical bias=0.057440, variance=0.010809, MSE=0.014108
  exact sampling bias=0.057143, variance=0.010714, MSE=0.013980
p=0.95, MLE: empirical bias=-0.000250, variance=0.004723, MSE=0.004723
  exact sampling bias=0.000000, variance=0.004750, MSE=0.004750
p=0.95, posterior mean: empirical bias=-0.128750, variance=0.002410, MSE=0.018986
  exact sampling bias=-0.128571, variance=0.002423, MSE=0.018954
A posterior mean can be assessed by repeated-sample risk; its Bayesian interpretation does not guarantee lower risk at every fixed p.

Compute both exact risks at both rates. The finite simulation estimates those risks; it cannot establish a universal ordering. Explain why the posterior mean is also a frequentist estimator rule.

Lab 3 · Numerical and modelling repairs

Download lab3_underflow_normalisation_and_prior.py

"""Diagnose probability underflow, parameter-dependent factors, and reused data."""
import math


def bernoulli_log_likelihood(p, k, n):
    if not 0 <= p <= 1 or not 0 <= k <= n:
        raise ValueError("invalid model or count")
    if (p == 0 and k > 0) or (p == 1 and n-k > 0):
        return -math.inf
    return (k*math.log(p) if k else 0.) + ((n-k)*math.log1p(-p) if n-k else 0.)


if __name__ == "__main__":
    n, k = 10000, 8000
    ll = [bernoulli_log_likelihood(p, k, n) for p in (.79, .8)]
    print(f"n={n}, k={k}: raw likelihoods={[math.exp(x) for x in ll]}")
    print(f"log likelihoods={ll[0]:.6f}, {ll[1]:.6f}; relative L(.79)/L(.8)={math.exp(ll[0]-ll[1]):.6f}")
    assert ll[1] > ll[0] and all(math.exp(x) == 0 for x in ll)
    for sigma in (1., 2.):
        full = -2*math.log(sigma)-math.log(2*math.pi)
        print(f"two zero residuals, sigma={sigma:.0f}: Gaussian full logL={full:.6f}; omitted-scale expression=0")
    a, b, k, n = 2, 2, 8, 10
    correct = (a+k)/(a+b+n)
    doubled = (a+2*k)/(a+b+2*n)
    print(f"One update: Beta(10,4), mean={correct:.6f}")
    print(f"Same data reused as new evidence: Beta(18,6), mean={doubled:.6f} (invalid independent-evidence claim)")
    sigma2, tau2, n, observed_mean = 1., 4., 4, 3.
    lam = sigma2/tau2
    map_mean = n*observed_mean/(n+lam)
    print(f"Gaussian mean, prior N(0,4), known noise variance1: MAP={map_mean:.6f}; sum-SSE penalty coefficient={lam:.6f}")
    print(f"For average-SSE convention the matching coefficient is lambda/n={lam/n:.6f}.")

n=10000, k=8000: raw likelihoods=[0.0, 0.0]
log likelihoods=-5007.074165, -5004.024235; relative L(.79)/L(.8)=0.047362
two zero residuals, sigma=1: Gaussian full logL=-1.837877; omitted-scale expression=0
two zero residuals, sigma=2: Gaussian full logL=-3.224171; omitted-scale expression=0
One update: Beta(10,4), mean=0.714286
Same data reused as new evidence: Beta(18,6), mean=0.750000 (invalid independent-evidence claim)
Gaussian mean, prior N(0,4), known noise variance1: MAP=2.823529; sum-SSE penalty coefficient=0.250000
For average-SSE convention the matching coefficient is lambda/n=0.062500.

Repair raw-product underflow, the omitted scale normalisation, reused evidence and the sum/average penalty coefficient. Name which parameter-dependent factors each optimisation needs.

Exercises 1–12 are required. Optional Exercises 13–14 add 35 minutes beyond the twelve-hour schedule.

Use MIT’s classical estimation lesson and Bayesian inference lesson. Mathematics for Machine Learning, Chapter 8, supplies the model-to-objective connection. The elementary estimators and beta update above are derived here.

Derive both MLEs, the variance denominator identity, MSE decomposition, beta normalisation/prediction and Gaussian-prior penalty scale. Explain identifiability, boundary existence and the three uncertainty objects without relying on a plot.

Exit task: compare MLE .8, posterior mean 5/7 and mode .75 for the ten-trial dataset. Explain their different objectives, then repair a duplicated update and a penalty whose coefficient was unchanged after averaging the residual loss.

Ready to move on: you can derive and interpret a point estimate under an explicit model. Inference next adds confidence coverage, hypothesis tests and experimental protocols. See the course overview.