# Estimation and Maximum Likelihood Concept Card: [Your study question or course] > A concept card for point estimation: estimators and estimates, bias, > variance, mean squared error, consistency and efficiency, the method > of moments, and maximum likelihood (likelihood, log-likelihood, score, > analytic MLEs for Bernoulli, Poisson and normal models, numerical MLE > with R `optim` and Python `scipy.optimize`, standard errors from the > observed information, the likelihood ratio test and AIC). It is for > students in a mathematical statistics or inference course, master's > and PhD students who fit models by maximum likelihood in a thesis, and > analysts who want to know what their software is doing when it says > "log-likelihood". It speeds up learning by pairing every formula with > runnable code on the same tiny data set, a fully worked example, two > faded ones, and a self-test with a separate answer key. The worked > examples use fictional data; replace them with your own. > Source basis: Fisher (1922, 1925), Pearson (1894), Wilks (1938), > Akaike (1974) and the NIST/SEMATECH e-Handbook, with all text and > examples written fresh. ## Document control | Field | Value | |---|---| | Course or project | [STAT 320 Mathematical Statistics / thesis methods chapter] | | Learner or analyst | [Your name] | | Data set | [Harbour City Library daily reports, n = 10 per branch; help-desk waiting times, n = 15 (fictional)] | | Software and versions | R 4.5.1 (stats); Python 3.14, numpy 2.5.4, scipy 1.18.1, statsmodels 0.15.0 | | Date | YYYY-MM-DD | | Version | [0.1] | ## The idea in one paragraph An estimator is a rule that turns data into a guess about an unknown parameter; the number it gives for your sample is the estimate. Good rules are, on average, close to the truth (low bias), do not jump around much from sample to sample (low variance), and home in on the truth as the sample grows (consistency). Maximum likelihood is the most widely used rule: choose the parameter value under which the observed data would have been most probable. It gives sensible estimates for almost any model you can write down, comes with standard errors from the curvature of the log-likelihood, and leads directly to likelihood ratio tests and AIC. Logistic, Poisson and mixed models are all fitted this way. ## Intuition Imagine 14 of 20 Lakemoor University (fictional) students pass a statistics quiz, and you want the pass probability p. Try candidate values: if p were 0.4, getting 14 passes would be fairly unlikely; if p were 0.95, getting as many as 6 failures would be very unlikely. The value that makes the observed result most likely is p = 0.7, the peak of the log-likelihood curve below. A sharp peak means the data rule out nearby values quickly (precise estimate, small SE); a flat peak means many values fit almost equally well (large SE). That curvature is the "information" in the data. A dartboard picture separates bias and variance. Bias is how far the centre of your cloud of darts is from the bullseye; variance is how spread out the cloud is. MSE combines the two, and sometimes a slightly biased thrower with a tight cloud beats an unbiased thrower with a wide one, as the simulation later shows. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Bernoulli log-likelihood for 14 passes in 20 trials, computed in R", "width": 380, "height": 220, "data": {"values": [ {"p": 0.40, "ll": -15.893}, {"p": 0.45, "ll": -14.766}, {"p": 0.50, "ll": -13.863}, {"p": 0.55, "ll": -13.161}, {"p": 0.60, "ll": -12.649}, {"p": 0.65, "ll": -12.330}, {"p": 0.70, "ll": -12.217}, {"p": 0.75, "ll": -12.345}, {"p": 0.80, "ll": -12.781}, {"p": 0.85, "ll": -13.658}, {"p": 0.90, "ll": -15.291}, {"p": 0.95, "ll": -18.692} ]}, "mark": {"type": "line", "point": true}, "encoding": { "x": {"field": "p", "type": "quantitative", "title": "Candidate pass probability p"}, "y": {"field": "ll", "type": "quantitative", "title": "Log-likelihood", "scale": {"zero": false}} } } ``` ```mermaid flowchart TD A["Choose a probability model for the data"] --> B["Write the likelihood L(theta)"] B --> C["Take logs: l(theta)"] C --> D{"Closed-form solution?"} D -- Yes --> E["Set the score to zero and solve"] D -- No --> F["Maximise numerically: optim or scipy.optimize"] E --> G["Check it is a maximum and inside the parameter space"] F --> G G --> H["SE from observed information (Hessian)"] H --> I["Wald CI, likelihood ratio test, AIC"] I --> J["Check model fit before trusting any of it"] ``` ## Formal definition ### Estimators and their properties An estimator $\hat{\theta}$ is a function of the random sample, so it has its own sampling distribution. Its bias is the gap between its average and the truth; its variance is its spread; MSE combines them. $$\text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta, \qquad \text{MSE}(\hat{\theta}) = E\left[(\hat{\theta} - \theta)^2\right] = \text{Var}(\hat{\theta}) + \text{Bias}(\hat{\theta})^2$$ - **Unbiased:** bias is zero for every θ. - **Consistent:** $\hat{\theta}$ gets arbitrarily close to θ with high probability as n grows. - **Efficient:** among unbiased estimators, it has the smallest variance. The Cramér-Rao bound says no unbiased estimator can have variance below $1/I(\theta)$, where $I(\theta)$ is the Fisher information. ### Method of moments Set sample moments equal to their population formulas and solve for the parameters (Pearson 1894). It is quick and often a good starting value for numerical MLE. For a gamma distribution with shape k and rate r, the mean is k/r and the variance is k/r², so $$\hat{k}_{MM} = \frac{\bar{x}^2}{\hat{\sigma}^2}, \qquad \hat{r}_{MM} = \frac{\bar{x}}{\hat{\sigma}^2}, \qquad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2$$ ### Maximum likelihood The likelihood is the joint probability (or density) of the observed data, read as a function of the parameter. Logs turn the product into a sum, which is easier to maximise and numerically stable. $$L(\theta) = \prod_{i=1}^n f(x_i \mid \theta), \qquad \ell(\theta) = \sum_{i=1}^n \log f(x_i \mid \theta), \qquad \hat{\theta}_{ML} = \arg\max_{\theta} \ell(\theta)$$ The score is the slope of the log-likelihood; the MLE usually solves score = 0. The observed information is the negative curvature at the peak; its inverse estimates the variance of the MLE. $$U(\theta) = \frac{\partial \ell}{\partial \theta}, \qquad J(\hat{\theta}) = -\frac{\partial^2 \ell}{\partial \theta^2}\Big|_{\hat{\theta}}, \qquad SE(\hat{\theta}) \approx \frac{1}{\sqrt{J(\hat{\theta})}}$$ With several parameters, $J$ is the Hessian matrix of $-\ell$ and the SEs are the square roots of the diagonal of its inverse. Under regular conditions (parameter not on a boundary, model correct, n reasonably large) the MLE is consistent, approximately normal and asymptotically efficient, and it is invariant: the MLE of $g(\theta)$ is $g(\hat{\theta})$. ### Analytic MLEs **Bernoulli** (x successes in n trials): $\ell(p) = x\log p + (n - x)\log(1-p)$. Setting $U(p) = x/p - (n-x)/(1-p) = 0$ gives $$\hat{p} = \frac{x}{n}, \qquad J(\hat{p}) = \frac{n}{\hat{p}(1-\hat{p})}, \qquad SE(\hat{p}) = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$ **Poisson** (counts $x_1, \dots, x_n$): $\ell(\lambda) = \sum x_i \log\lambda - n\lambda - \sum \log x_i!$. Setting $U(\lambda) = \sum x_i/\lambda - n = 0$ gives $$\hat{\lambda} = \bar{x}, \qquad J(\hat{\lambda}) = \frac{\sum x_i}{\hat{\lambda}^2} = \frac{n}{\hat{\lambda}}, \qquad SE(\hat{\lambda}) = \sqrt{\frac{\hat{\lambda}}{n}}$$ **Normal** (mean μ, variance σ²): solving both score equations gives $$\hat{\mu} = \bar{x}, \qquad \hat{\sigma}^2_{ML} = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2$$ The MLE of σ² divides by n, not n − 1, so it is biased downwards by the factor (n − 1)/n. It is still consistent, because that factor tends to 1. For the 12 Brightwell Bakery loaf weights used in other cards, R gives $\hat{\mu}$ = 497.55, $\hat{\sigma}_{ML}$ = 2.1654 and s = 2.2617. ### Likelihood ratio test and AIC To compare a simpler (nested) model with a fuller one, compare their maximised log-likelihoods (Wilks 1938). Under H0 (the simpler model is true) the statistic is approximately chi-square with df equal to the number of extra parameters. $$\Lambda = 2\left[\ell(\hat{\theta}_{full}) - \ell(\hat{\theta}_{0})\right] \sim \chi^2_{df}$$ AIC (Akaike 1974) trades fit against complexity for any set of models fitted to the same data, nested or not; smaller is better. $$\text{AIC} = -2\,\ell(\hat{\theta}) + 2k$$ | Symbol | Meaning | |---|---| | $\theta$, $\hat{\theta}$ | Parameter and its estimator or estimate | | $L(\theta)$, $\ell(\theta)$ | Likelihood and log-likelihood | | $U(\theta)$ | Score: first derivative of $\ell$ | | $J(\hat{\theta})$ | Observed information: minus the second derivative at the MLE | | $I(\theta)$ | Expected (Fisher) information | | $p$, $\lambda$, $\mu$, $\sigma^2$ | Bernoulli, Poisson and normal parameters | | $k$, $r$ | Gamma shape and rate | | $\Lambda$ | Likelihood ratio statistic | | $k$ in AIC | Number of estimated parameters | | $n$ | Sample size n | ## Worked example A: fully solved **Question.** Harbour City Library (fictional) logs lost-property reports per day. Over ten days the Central branch had the counts below. Estimate the daily rate, give an SE and 95% CI, and test whether the Harbour branch has the same rate. | Day | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |---|---|---|---|---|---|---|---|---|---|---| | Central | 3 | 5 | 2 | 4 | 6 | 3 | 4 | 5 | 2 | 4 | | Harbour | 1 | 2 | 0 | 3 | 1 | 2 | 2 | 1 | 0 | 2 | **Step 1: model.** Counts of independent events per day at a constant rate: Poisson. **Step 2: analytic MLE.** From the derivation, $\hat{\lambda} = \bar{x}$ = 38 / 10 = 3.8 reports per day. **Step 3: SE.** $SE = \sqrt{3.8/10}$ = 0.616. **Step 4: Wald CI.** 3.8 ± 1.96 × 0.616 = [2.59, 5.01]. **Step 5: likelihood ratio test.** H0: both branches share one rate (one parameter, $\hat{\lambda}$ = 52 / 20 = 2.6). Full model: separate rates (two parameters, 3.8 and 1.4). One extra parameter, so df = 1. Outputs on this card are trimmed to the key lines (for example the header of R's analysis-of-deviance table is dropped). ```r central <- c(3, 5, 2, 4, 6, 3, 4, 5, 2, 4) harbour <- c(1, 2, 0, 3, 1, 2, 2, 1, 0, 2) lam <- mean(central); c(lambda_hat = lam, se = sqrt(lam / length(central))) ll0 <- sum(dpois(c(central, harbour), mean(c(central, harbour)), log = TRUE)) ll1 <- sum(dpois(central, mean(central), log = TRUE)) + sum(dpois(harbour, mean(harbour), log = TRUE)) LR <- 2 * (ll1 - ll0) round(c(ll0 = ll0, ll1 = ll1, LR = LR, p = pchisq(LR, 1, lower.tail = FALSE)), 5) # The same test as a Poisson GLM branch <- rep(c("central", "harbour"), each = 10) anova(glm(c(central, harbour) ~ branch, family = poisson), test = "LRT") ``` ```text lambda_hat se 3.8000000 0.6164414 ll0 ll1 LR p -37.53596 -31.78191 11.50811 0.00069 Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 19 23.998 branch 1 11.508 18 12.490 0.0006929 *** ``` ```python import numpy as np import statsmodels.api as sm from scipy import stats central = np.array([3, 5, 2, 4, 6, 3, 4, 5, 2, 4]) harbour = np.array([1, 2, 0, 3, 1, 2, 2, 1, 0, 2]) print(central.mean(), round(np.sqrt(central.mean() / 10), 4)) y = np.concatenate([central, harbour]); grp = np.repeat([0, 1], 10) m0 = sm.GLM(y, np.ones(20), family=sm.families.Poisson()).fit() m1 = sm.GLM(y, sm.add_constant(grp), family=sm.families.Poisson()).fit() LR = 2 * (m1.llf - m0.llf) print(round(m0.llf, 5), round(m1.llf, 5), round(LR, 5), round(stats.chi2.sf(LR, 1), 5)) ``` ```text 3.8 0.6164 -37.53596 -31.78191 11.50811 0.00069 ``` **Step 6: interpret.** The Central branch logs about 3.8 reports a day (95% CI 2.59 to 5.01). Allowing separate rates raises the log-likelihood by 5.75, giving Λ = 11.51 on 1 df, p < .001: the branches' rates differ. The GLM deviance drop (23.998 − 12.490 = 11.508) is exactly the same LR statistic, because deviance is −2 times the log-likelihood up to a constant. *Report:* A Poisson model estimated 3.8 lost-property reports per day at Central (95% CI [2.59, 5.01]) and 1.4 at Harbour; a likelihood ratio test favoured separate rates, χ²(1) = 11.51, *p* < .001. **Check before trusting it.** The Poisson model assumes the variance equals the mean. If counts are more spread out (overdispersion), the SEs are too small; see stats-count-regression. ## Worked example B: half solved **Question.** Fifteen waiting times (minutes) at Lakemoor University's IT help desk (fictional): 2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3. Waiting times are positive and right-skewed, so fit a gamma distribution. The gamma shape has no closed-form MLE, so we maximise numerically, starting from the method of moments. ```r x <- c(2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3) v <- mean((x - mean(x))^2) mom <- c(shape = mean(x)^2 / v, rate = mean(x) / v); round(mom, 4) # Work on the log scale so that optim cannot propose negative values nll <- function(par) -sum(dgamma(x, shape = exp(par[1]), rate = exp(par[2]), log = TRUE)) fit <- optim(log(mom), nll, method = "BFGS", hessian = TRUE) est <- exp(fit$par); round(c(est, logLik = -fit$value), 4) se <- est * sqrt(diag(solve(fit$hessian))) # delta method back to the natural scale round(se, 4) ll_exp <- sum(dexp(x, 1 / mean(x), log = TRUE)) # exponential = gamma with shape 1 round(c(AIC_exp = -2 * ll_exp + 2, AIC_gamma = 2 * fit$value + 4), 3) LR <- 2 * (-fit$value - ll_exp); round(c(LR = LR, p = pchisq(LR, 1, lower.tail = FALSE)), 4) ``` ```text shape rate 3.1734 0.8054 shape rate logLik 3.4524 0.8762 -31.0026 shape rate 1.2047 0.3291 AIC_exp AIC_gamma 73.135 66.005 LR p 9.1302 0.0025 ``` ```python import numpy as np from scipy import stats, optimize from statsmodels.tools.numdiff import approx_hess x = np.array([2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3]) v = x.var() # ddof=0, the moment estimator mom = np.array([x.mean()**2 / v, x.mean() / v]); print(np.round(mom, 4)) def nll(par): shape, rate = par if shape <= 0 or rate <= 0: return np.inf return -stats.gamma.logpdf(x, a=shape, scale=1 / rate).sum() fit = optimize.minimize(nll, mom, method="Nelder-Mead", options={"xatol": 1e-8, "fatol": 1e-10}) print(np.round(fit.x, 4), round(-fit.fun, 4)) H = approx_hess(fit.x, nll) # observed information matrix print(np.round(np.sqrt(np.diag(np.linalg.inv(H))), 4)) a, loc, scale = stats.gamma.fit(x, floc=0); print(round(a, 4), round(1 / scale, 4)) ll_exp = stats.expon.logpdf(x, scale=x.mean()).sum() print("AIC exp", round(-2 * ll_exp + 2, 3), "AIC gamma", round(2 * fit.fun + 4, 3)) ``` ```text [3.1734 0.8054] [3.4524 0.8762] -31.0026 [1.2047 0.3291] 3.4524 0.8762 AIC exp 73.135 AIC gamma 66.005 ``` Note the parameterisations: R's `dgamma` takes a rate, while scipy's `gamma` takes a scale = 1 / rate; scipy's `fit` also estimates a location unless you fix `floc=0`. Mixing them up is a classic error. 1. Method-of-moments estimates: shape = [ ], rate = [ ]. 2. MLE: shape = [ ], rate = [ ]; maximised log-likelihood = [ ]. 3. Implied mean waiting time = shape / rate = [ ] minutes. Compare it with $\bar{x}$ = 3.94: [same / different]? (Hint: the gamma MLE always matches the sample mean.) 4. 95% Wald CI for the shape: 3.4524 ± 1.96 × [ ] = [ ]. 5. Exponential versus gamma: Λ = [ ] on [ ] df, p = [ ]; AIC prefers [ ] by [ ] points. 6. Report: "A gamma model (shape [ ], rate [ ] per minute) fitted the waiting times better than an exponential model, χ²(1) = [ ], *p* = [ ], AIC [ ] versus [ ]." ## Exercise C: unsolved 1. Brightwell Bakery (fictional) records the minutes between customer arrivals: 2, 4, 3, 7, 4. Assume an exponential model with density $f(x) = \lambda e^{-\lambda x}$. Derive the MLE of λ, compute it, and give its SE from the observed information. 2. Compute the maximised log-likelihood. 3. For the quiz example (14 passes in 20), compute the 95% Wald CI for p and the likelihood ratio statistic for H0: p = 0.5 (use the log-likelihood values −12.217 at 0.7 and −13.863 at 0.5). 4. By the invariance property, what is the MLE of the odds p / (1 − p)? ## Simulation you can run **Bias versus MSE.** Draw 20,000 samples of n = 5 from a normal population with σ² = 4. Compare the MLE of the variance (divide by n) with the unbiased version (divide by n − 1). ```r set.seed(2026) sims <- replicate(20000, { y <- rnorm(5, 0, 2) c(mle = mean((y - mean(y))^2), unbiased = var(y)) }) round(rowMeans(sims), 3) # average estimate (truth = 4) round(rowMeans((sims - 4)^2), 3) # MSE ``` ```text mle unbiased 3.208 4.010 mle unbiased 5.750 8.004 ``` ```python import numpy as np rng = np.random.default_rng(2026) ys = rng.normal(0, 2, size=(20000, 5)) mle, unb = ys.var(axis=1, ddof=0), ys.var(axis=1, ddof=1) print(np.round([mle.mean(), unb.mean()], 3), np.round([((mle - 4)**2).mean(), ((unb - 4)**2).mean()], 3)) ``` ```text [3.188 3.986] [5.724 7.915] ``` Theory predicts an average of $\frac{n-1}{n}\sigma^2$ = 3.2 for the MLE and 4 for the unbiased estimator, with MSEs of 5.76 and 8.00. Both simulations land close (R and Python use different random generators, so their numbers differ slightly). The lesson: the MLE is biased, yet it has the smaller mean squared error, because its lower variance more than pays for its bias. "Unbiased" is not the same as "best". ## Misconceptions | Wrong statement | Correct statement | |---|---| | The likelihood is the probability that the parameter equals θ. | The likelihood is the probability of the observed data given θ, viewed as a function of θ. It does not integrate to 1 over θ and is not a distribution for θ (that needs a prior; see stats-bayesian-inference). | | Unbiased estimators are always better. | MSE matters; a small bias with much lower variance can be preferable (the n versus n − 1 simulation). | | MLEs are unbiased. | Many are biased in small samples (σ² divides by n); they are consistent and asymptotically unbiased. | | A higher log-likelihood always means a better model. | Adding parameters never lowers the maximised log-likelihood; use a likelihood ratio test for nested models or AIC to penalise complexity. | | AIC values can be compared across different data sets. | AIC is only comparable between models fitted to exactly the same observations and response. | | The Hessian-based SE is exact. | It is a large-sample approximation; for small n or parameters near a boundary, prefer profile-likelihood or bootstrap intervals. | | If optim stops, it found the maximum. | Check the convergence code, try different starting values, and confirm the score is near zero; likelihoods can have several peaks or flat ridges. | | The likelihood ratio test works for any comparison. | It needs nested models and parameters away from the boundary of the parameter space (testing a variance = 0 breaks the chi-square reference). | ## Connections to other concepts - **Distributions** (stats-distributions-sheet): every likelihood starts from a distribution's pmf or pdf. - **Sampling distributions** (stats-sampling-distributions): an estimator's bias and variance describe its sampling distribution. - **Confidence intervals** (stats-confidence-intervals): Wald intervals are estimate ± 1.96 × SE from the information; profile-likelihood intervals invert the LR test. - **Hypothesis testing** (stats-hypothesis-testing): Wald, score and likelihood ratio tests are three large-sample tests built from the same log-likelihood. - **Regression models** (stats-logistic-regression, stats-count-regression, stats-mixed-models): all fitted by maximum likelihood (or REML), and their output reports log-likelihood, AIC and deviance. - **Bayesian inference** (stats-bayesian-inference): posterior ∝ prior × likelihood; with a flat prior the posterior mode equals the MLE. ## Self-test 1. Define bias, variance and MSE, and state the identity that links them. 2. Derive the MLE of p for x successes in n Bernoulli trials. 3. Why does a sharper log-likelihood peak imply a smaller SE? 4. What is the MLE of σ² for a normal sample, is it biased, and is it consistent? 5. Two nested models have log-likelihoods −120.4 (3 parameters) and −116.1 (5 parameters). Compute Λ, its df, and the AIC of each model. 6. Give one reason to start a numerical optimiser from method-of-moments estimates. ## My notes and next steps | Item | Entry | |---|---| | Confidence today (1-5) | [3] | | Hardest idea for me | [*Example:* observed versus expected information] | | Open questions | [Question for the tutor or supervisor] | | Mistakes logged in stats-error-log | [Entry IDs] | | Next review date (spacing) | YYYY-MM-DD | | Next template to study | [stats-bayesian-inference] | Derivation I can now do from memory: [list them, for example Poisson MLE and its SE]. ## Sources and further reading Structure aligned with a standard first course in mathematical statistics and with the GAISE College Report (ASA 2016) emphasis on simulation; based on the original papers below. All wording and examples are original. - Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19(6), 716–723. - Fisher, R. A. (1922). On the mathematical foundations of theoretical statistics. *Philosophical Transactions of the Royal Society of London, Series A*, 222, 309–368. - Fisher, R. A. (1925). Theory of statistical estimation. *Proceedings of the Cambridge Philosophical Society*, 22(5), 700–725. - Pearson, K. (1894). Contributions to the mathematical theory of evolution. *Philosophical Transactions of the Royal Society of London, Series A*, 185, 71–110. - Wilks, S. S. (1938). The large-sample distribution of the likelihood ratio for testing composite hypotheses. *Annals of Mathematical Statistics*, 9(1), 60–62. - Nelder, J. A., and Mead, R. (1965). A simplex method for function minimization. *The Computer Journal*, 7(4), 308–313. - NIST/SEMATECH e-Handbook of Statistical Methods, section 1.3.6.5.2 Maximum Likelihood, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain. - GAISE College Report ASA Revision Committee (2016). *Guidelines for Assessment and Instruction in Statistics Education (GAISE) College Report 2016*. American Statistical Association. - Next steps: Pawitan, Y. (2001). *In All Likelihood: Statistical Modelling and Inference Using Likelihood*. Oxford University Press; Casella, G., and Berger, R. L. (2002). *Statistical Inference* (2nd ed.). Duxbury. - Further reading: de Smith, M. J. (2024). *Statistical Analysis Handbook*, 2024 edition. StatsRef, https://www.statsref.com/; OpenIntro Statistics (CC BY-SA). ## Answer key **Example B.** 1. Shape 3.1734, rate 0.8054. 2. Shape 3.4524, rate 0.8762; log-likelihood −31.0026. 3. 3.4524 / 0.8762 = 3.94 minutes, the same as $\bar{x}$. 4. 3.4524 ± 1.96 × 1.2047 = [1.09, 5.81]. 5. Λ = 9.13 on 1 df, p = .0025; AIC prefers the gamma model by 7.13 points (66.005 versus 73.135). 6. "A gamma model (shape 3.45, rate 0.88 per minute) fitted the waiting times better than an exponential model, χ²(1) = 9.13, *p* = .003, AIC 66.0 versus 73.1." **Exercise C.** 1. $\ell(\lambda) = n\log\lambda - \lambda\sum x_i$; $U = n/\lambda - \sum x_i = 0$ gives $\hat{\lambda} = 1/\bar{x}$ = 1 / 4 = 0.25 per minute. $J = n/\hat{\lambda}^2$, so $SE = \hat{\lambda}/\sqrt{n}$ = 0.25 / √5 = 0.112. 2. $5\log 0.25 - 0.25 \times 20$ = −11.93. 3. Wald CI 0.7 ± 1.96 × 0.1025 = [0.499, 0.901]. Λ = 2 × (−12.217 − (−13.863)) = 3.29 on 1 df, p = 0.070: not significant at α = 0.05. 4. 0.7 / 0.3 = 2.33. **Self-test.** 1. Bias = E[estimator] − truth; variance = spread of the estimator across samples; MSE = average squared error = variance + bias². 2. $\ell(p) = x\log p + (n-x)\log(1-p)$; derivative $x/p - (n-x)/(1-p) = 0$ gives $\hat{p} = x/n$. 3. The SE is roughly one over the square root of the curvature; a sharp peak means high curvature (much information), so the SE is small. 4. $\frac{1}{n}\sum(x_i - \bar{x})^2$; biased (expectation $\frac{n-1}{n}\sigma^2$) but consistent. 5. Λ = 2 × (−116.1 + 120.4) = 8.6 on 2 df (p ≈ 0.014). AIC = 246.8 (3 parameters) and 242.2 (5 parameters): the larger model is preferred. 6. Good starting values speed convergence and reduce the risk of the optimiser wandering into an invalid region or a local optimum. ## Revision history | Version | Date | Author | Change summary | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Your name] | First draft from template | | 0.2 | YYYY-MM-DD | [Your name] | *Example:* Added Weibull fit for my thesis data |