HW #5: Modeling an Unknown Distribution
Dr. Teresa Vasquez, EdD, MBA, MSDS (in progress)
October 01, 2026

This report was compiled from my own Jupyter notebook (HW5.ipynb, included in this folder) with a custom Claude Code script built by Dr. Teresa Vasquez. The math, code, and answers are my worked solutions; the script converts the notebook into this R Markdown file and knits it. Every number below is recomputed by R when the file is knit.

Suppose you have 1,000 measurements and want a formula that describes where they came from. That formula is a probability distribution, and it has dials, called parameters, that need to be set. This tutorial walks through two common ways to set them, the method of moments (MM) and maximum likelihood estimation (MLE), using the heights of 1,000 adult women from the National Health and Nutrition Examination Survey (NHANES) 2009 to 2010. We try three candidate distributions (Normal, Gamma, and Exponential), check each fit with plots, estimate the median, and measure how much that median could vary from sample to sample.

Vocabulary

Term What it means here Textbook (Matloff)
Parameter A dial that sets the shape of a distribution, like the center or the spread Section 5.2, page 86
Probability density function (PDF) The smooth curve showing where values are likely to fall Section 6.5, page 119
Histogram Bars showing how many heights fall in each range Section 8.2, page 172
Cumulative distribution function (CDF) For any height, the fraction of women at or below it Section 6.4.1, page 115
Empirical cumulative distribution function (ECDF) The same idea, counted straight from the data Section 8.6, page 187
Normal distribution The bell curve, set by its center and spread Section 6.7.2, page 127
Exponential distribution A curve that starts high at 0 and falls, built for waiting times Section 6.7.3, page 128
Gamma distribution A curve for positive values that can lean to one side Section 6.7.4, page 131
Method of moments (MM) Set the model’s averages equal to the data’s averages and solve Section 8.4.1, page 181
Likelihood A score for how probable a choice of dials makes the data we saw Section 8.4.3, page 183
Maximum likelihood estimation (MLE) Pick the dials with the highest likelihood score Section 8.4.3, page 183
Solving numerically Letting the computer find an answer by trial and improvement when no formula exists Section 8.4.3, page 183
Unbiased Right on average; why the variance divides by n − 1 Section 7.6.1, page 159
Kolmogorov-Smirnov (KS) distance The biggest vertical gap between the model’s cumulative curve and the data’s Section 8.6, page 187
Median The height that splits the women in half
Sampling distribution How an estimate would vary if we drew a new sample Section 7.2, page 150
Standard error The typical size of that variation Section 7.7, page 161
Simulation Using the computer to generate many fake samples Chapter 2, page 35
Quantile-quantile (QQ) plot Sorted data plotted against where the model says each value should fall; points on the diagonal mean a match Named in the assignment

Get Started

The three distributions and how R names their dials:

Family Parameters Support Mean Variance R functions
Normal \(\mu, \sigma\) all real numbers \(\mu\) \(\sigma^2\) dnorm pnorm qnorm rnorm
Gamma shape \(\alpha\), rate \(\beta\) \(x > 0\) \(\alpha/\beta\) \(\alpha/\beta^2\) dgamma pgamma qgamma rgamma
Exponential rate \(\lambda\) \(x > 0\) \(1/\lambda\) \(1/\lambda^2\) dexp pexp qexp rexp 
# Setup: seed, packages, plot size, and the black-background plot style
set.seed(6110)
suppressPackageStartupMessages(library(tidyverse))
options(repr.plot.width = 15, repr.plot.height = 9)
dark <- function() par(bg = "black", col.axis = "white", col.lab = "white", col.main = "white", fg = "white")
ink <- "white"   # color for the empirical cumulative distribution function steps

Part 0: Load the data and get your raw ingredients

Before fitting anything, compute a few summary numbers. Every estimate later in this tutorial is built from them.

\(x\) = the height in centimeters of each of 1,000 adult women.

Step 1. Load the data

d1.csv was written from the homework’s exact code (shown in the code cell). length(x) \(= 1000\) and sum(is.na(x)) \(= 0\), so there are 1,000 heights and no missing values. Units are centimeters.

# d1.csv was created once with the homework's code:
#   Hmisc::getHdata(nhgh)
#   d1 <- nhgh %>% filter(sex == "female") %>% filter(age >= 18) %>%
#     select(ht) %>% filter(1:n() <= 1000)
#   write.csv(d1, "d1.csv", row.names = FALSE)
if (file.exists("d1.csv")) {
  d1 <- read.csv("d1.csv")
} else {
  Hmisc::getHdata(nhgh)
  d1 <- nhgh %>% filter(sex == "female") %>% filter(age >= 18) %>%
    select(ht) %>% filter(1:n() <= 1000)
}
x <- d1$ht
n <- length(x)
c(n = n, missing = sum(is.na(x)))
##       n missing 
##    1000       0

Step 2. Find the average height

\[\bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i = \frac{160741.9}{1000} = 160.7419 \text{ cm}\]

xbar <- mean(x)
c(sum = sum(x), xbar = xbar)
##         sum        xbar 
## 160741.9000    160.7419

Step 3. Measure the spread two ways

Divide by \(n\):

\[\hat\sigma^2_{n} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2 = \frac{53531.1744}{1000} = 53.5312 \text{ cm}^2, \qquad \sqrt{53.5312} = 7.3165 \text{ cm}\]

Divide by \(n - 1\) (what var() does):

\[s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2 = \frac{53531.1744}{999} = 53.5848 \text{ cm}^2, \qquad s = \sqrt{53.5848} = 7.3202 \text{ cm}\]

The two differ by the factor \(1000/999 = 1.001\), about 0.1 percent.

Choice: divide by \(n - 1\) (var(x) \(= 53.5848\)) for the method of moments (MM).

This is the standard sample variance, the one R’s var() computes. The \(n - 1\) version is unbiased for \(\sigma^2\).

The \(1/n\) version \(53.5312\) is not used for the method of moments (MM); it shows up on its own as the Normal maximum likelihood estimation (MLE) in Part 1, which is biased by the factor \((n-1)/n\).

ss <- sum((x - xbar)^2)
s2 <- mean((x - xbar)^2)      # divide by n: what the Normal maximum likelihood estimation (MLE) produces
v <- var(x)                   # divide by n - 1: used for method of moments (MM), the standard sample variance
c(sum_sq = ss, s2_n = s2, var_n_minus_1 = v, s_n = sqrt(s2), sd_n_minus_1 = sqrt(v),
  ratio = v / s2)
##        sum_sq          s2_n var_n_minus_1           s_n  sd_n_minus_1 
##  53531.174390     53.531174     53.584759      7.316500      7.320161 
##         ratio 
##      1.001001

Step 4. Average the log heights

\[\overline{\log x} = \frac{1}{n}\sum_{i=1}^{n}\log x_i = 5.078763\]

For comparison, \(\log \bar{x} = \log(160.7419) = 5.079800\). The gap \(\log\bar{x} - \overline{\log x} = 0.001037\) is tiny because the heights are tightly clustered relative to their size. This gap drives the Gamma shape estimate in Part 3.

mlog <- mean(log(x))
c(mean_log_x = mlog, log_xbar = log(xbar), gap = log(xbar) - mlog)
##  mean_log_x    log_xbar         gap 
## 5.078763038 5.079799973 0.001036935

Step 5. Check the shape

\[\operatorname{median}(x) = 160.8, \qquad \min = 139.4, \qquad \max = 190.8 \text{ cm}\]

Mean minus median \(= 160.7419 - 160.8 = -0.0581\) cm, a tiny gap compared with \(s = 7.3202\) cm. The distance from the mean to the minimum is \(160.7419 - 139.4 = 21.34\) cm and to the maximum is \(190.8 - 160.7419 = 30.06\) cm, so the upper tail reaches a little farther. The sample skewness \(\frac{1}{n}\sum (x_i - \bar{x})^3 / \hat\sigma_n^3 = 0.0757\) is close to 0. Conclusion: the heights are close to symmetric, with at most a very slight right lean.

c(median = median(x), min = min(x), max = max(x),
  mean_minus_median = xbar - median(x),
  skewness = mean((x - xbar)^3) / s2^1.5)
##            median               min               max mean_minus_median 
##      160.80000000      139.40000000      190.80000000       -0.05810000 
##          skewness 
##        0.07566446

Prediction

Prediction: the Exponential will fit worst by far. Its density is highest at 0 cm and decreases from there, and it has one parameter that forces the standard deviation to equal the mean (about 161 cm), while the data are a symmetric hump near 161 cm with standard deviation about 7.3 cm. The Normal should fit well because the data look symmetric. The Gamma should look almost identical to the Normal, since a Gamma with a large shape is nearly bell shaped.

Part 1: Normal, by method of moments (MM) and maximum likelihood estimation (MLE)

The Normal distribution is the classic bell curve. It has two dials: μ, the center, and σ, the spread. We set them both ways and compare.

\(\mu\) is the typical woman’s height, the center of the bell. \(\sigma\) is how far a typical woman sits from that center.

Step 1. Method of moments (MM): match the averages

Model moments: \(E[X] = \mu\) and \(\operatorname{Var}(X) = \sigma^2\).

Moment equations:

\[\mu = \bar{x}, \qquad \sigma^2 = s^2 = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})^2\]

Solve:

\[\hat{\mu}_{\text{MM}} = \bar{x} = 160.7419 \text{ cm}\]

\[\hat{\sigma}_{\text{MM}} = \sqrt{s^2} = \sqrt{\frac{53531.1744}{999}} = \sqrt{53.5848} = 7.3202 \text{ cm}\]

(Subscript on the estimate: method of moments (MM).)

mu_mm <- xbar
sigma_mm <- sqrt(v)
c(mu_mm = mu_mm, sigma_mm = sigma_mm)
##      mu_mm   sigma_mm 
## 160.741900   7.320161

Step 2. Maximum likelihood estimation (MLE): write the log-likelihood

The likelihood scores a choice of dials by how probable it makes the data we actually saw. Taking the log turns a product of 1,000 tiny numbers into a sum, which is easier to work with and peaks at the same place.

One point:

\[\log f(x_i; \mu, \sigma) = \log\left(\frac{1}{\sigma\sqrt{2\pi}}\right) - \frac{(x_i - \mu)^2}{2\sigma^2} = -\log\sigma - \tfrac{1}{2}\log(2\pi) - \frac{(x_i - \mu)^2}{2\sigma^2}\]

Sum over \(n\) points:

\[\ell(\mu, \sigma) = \sum_{i=1}^{n}\log f(x_i; \mu, \sigma) = -n\log\sigma - \frac{n}{2}\log(2\pi) - \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i - \mu)^2\]

loglik_norm <- function(mu, sigma) -n * log(sigma) - n / 2 * log(2 * pi) - sum((x - mu)^2) / (2 * sigma^2)
c(formula = loglik_norm(160, 7), dnorm_sum = sum(dnorm(x, 160, 7, log = TRUE)))
##   formula dnorm_sum 
## -3416.702 -3416.702

Step 3. Solve for the center

Only the last term of \(\ell\) contains \(\mu\). By the chain rule, \(\frac{\partial}{\partial\mu}(x_i - \mu)^2 = -2(x_i - \mu)\):

\[\frac{\partial \ell}{\partial \mu} = -\frac{1}{2\sigma^2}\sum_{i=1}^{n} -2(x_i - \mu) = \frac{1}{\sigma^2}\sum_{i=1}^{n}(x_i - \mu)\]

Set to zero (\(\sigma^2 > 0\) so it can be multiplied away):

\[\sum_{i=1}^{n}(x_i - \mu) = 0 \;\Longrightarrow\; \sum_{i=1}^{n} x_i - n\mu = 0 \;\Longrightarrow\; \hat{\mu}_{\text{MLE}} = \frac{1}{n}\sum_{i=1}^{n} x_i = \bar{x}\]

Second derivative: \(\frac{\partial^2 \ell}{\partial \mu^2} = -\frac{n}{\sigma^2} < 0\), so this is a maximum.

\[\hat{\mu}_{\text{MLE}} = 160.7419 \text{ cm}\]

(Subscript on the estimate: maximum likelihood estimation (MLE).)

Step 4. Solve for the spread

\[\frac{\partial \ell}{\partial \sigma} = -\frac{n}{\sigma} - \frac{1}{2}\sum_{i=1}^{n}(x_i - \mu)^2 \cdot \frac{d}{d\sigma}\sigma^{-2} = -\frac{n}{\sigma} + \frac{1}{\sigma^3}\sum_{i=1}^{n}(x_i - \mu)^2\]

Set to zero and multiply both sides by \(\sigma^3\):

\[-n\sigma^2 + \sum_{i=1}^{n}(x_i - \mu)^2 = 0 \;\Longrightarrow\; \sigma^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu)^2\]

Plug in \(\hat{\mu} = \bar{x}\):

\[\hat{\sigma}^2_{\text{MLE}} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2 = \frac{53531.1744}{1000} = 53.5312\]

\[\hat{\sigma}_{\text{MLE}} = \sqrt{53.5312} = 7.3165 \text{ cm}\]

Second derivative check at the solution: \(\frac{\partial^2 \ell}{\partial \sigma^2} = \frac{n}{\sigma^2} - \frac{3}{\sigma^4}\sum(x_i - \bar{x})^2 = \frac{n}{\hat\sigma^2} - \frac{3n\hat\sigma^2}{\hat\sigma^4} = -\frac{2n}{\hat{\sigma}^2} < 0\), a maximum.

The maximized log-likelihood is \(\ell(\hat\mu, \hat\sigma) = -3409.0706\).

mu_mle <- xbar
sigma_mle <- sqrt(mean((x - mu_mle)^2))
num <- optim(c(150, 5), function(p) -loglik_norm(p[1], p[2]))
c(mu_mle = mu_mle, sigma_mle = sigma_mle, optim_mu = num$par[1], optim_sigma = num$par[2],
  loglik = loglik_norm(mu_mle, sigma_mle))
##       mu_mle    sigma_mle     optim_mu  optim_sigma       loglik 
##   160.741900     7.316500   160.739829     7.316427 -3409.070625

How the two methods compare

The centers are identical: \(\hat\mu = 160.7419\) from both. The spreads differ slightly: \(\hat\sigma_{\text{MM}} = 7.3202\) vs \(\hat\sigma_{\text{MLE}} = 7.3165\), a difference of \(0.0037\) cm (0.05 percent). The cause is only the divisor. The maximum likelihood estimation (MLE) calculus produces \(1/n\); the method of moments (MM) here uses \(n - 1\). Had the method of moments (MM) used \(1/n\), the two would match exactly. With \(n = 1000\) the difference is far too small to see in any plot.

Part 2: Gamma, by method of moments (MM)

The Gamma forbids negative heights and allows a lean to one side. Shape \(\alpha\) controls how bell-like it is; rate \(\beta\) controls the scale. Mean \(= \alpha/\beta\), variance \(= \alpha/\beta^2\).

Step 1. Write the moment equations

\[\text{(1)}\quad \frac{\alpha}{\beta} = \bar{x} = 160.7419\]

\[\text{(2)}\quad \frac{\alpha}{\beta^2} = s^2 = 53.5848 \quad (\text{divide by } n - 1)\]

Step 2. Solve for the rate

\[\frac{(1)}{(2)}: \quad \frac{\alpha/\beta}{\alpha/\beta^2} = \frac{\alpha}{\beta}\cdot\frac{\beta^2}{\alpha} = \beta = \frac{\bar{x}}{s^2}\]

\(\alpha\) cancels, and one factor of \(\beta\) cancels, leaving \(\beta\) alone.

\[\hat{\beta}_{\text{MM}} = \frac{\bar{x}}{s^2}\]

Step 3. Solve for the shape

From (1), \(\alpha = \beta\,\bar{x}\). Substitute \(\hat\beta = \bar{x}/s^2\):

\[\hat{\alpha}_{\text{MM}} = \frac{\bar{x}}{s^2}\cdot\bar{x} = \frac{\bar{x}^2}{s^2}, \qquad \hat{\beta}_{\text{MM}} = \frac{\bar{x}}{s^2}\]

Step 4. Plug in the numbers

\[\hat{\beta}_{\text{MM}} = \frac{160.7419}{53.5848} = 2.9998\]

\[\hat{\alpha}_{\text{MM}} = \frac{160.7419^2}{53.5848} = \frac{25837.96}{53.5848} = 482.19\]

Check: \(\hat\alpha/\hat\beta = 482.19 / 2.9998 = 160.7419\) (matches \(\bar{x}\)) and \(\hat\alpha/\hat\beta^2 = 53.5848\) (matches \(s^2\)).

(With the \(1/n\) variance instead, you would get \(\hat\alpha = 482.67\) and \(\hat\beta = 3.0028\).)

b_mm <- xbar / v
a_mm <- xbar^2 / v
c(a_mm = a_mm, b_mm = b_mm, model_mean = a_mm / b_mm, model_var = a_mm / b_mm^2)
##       a_mm       b_mm model_mean  model_var 
## 482.188571   2.999769 160.741900  53.584759

What a large shape means

A Gamma with shape \(\alpha\) has skewness \(2/\sqrt{\alpha}\). With \(\hat\alpha = 482.19\) that is \(2/\sqrt{482.19} = 0.0911\), almost 0, so this Gamma is a nearly symmetric bell. Its plots should be almost indistinguishable from the Normal’s, with the Gamma curve very slightly shifted so its right tail is a touch longer.

Part 3: Gamma, by maximum likelihood estimation (MLE)

Same Gamma family, different tuning rule: which \(\alpha\) and \(\beta\) make these exact 1,000 women the most plausible?

Step 1. Write the log-likelihood

Gamma density with shape \(\alpha\) and rate \(\beta\): \(f(x) = \dfrac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}\). Taking the log of one point:

\[\log f(x_i; \alpha, \beta) = \alpha\log\beta - \log\Gamma(\alpha) + (\alpha - 1)\log x_i - \beta x_i\]

Sum over \(n\) points, using \(\sum x_i = n\bar{x}\) and \(\sum \log x_i = n\,\overline{\log x}\):

\[\ell(\alpha, \beta) = n\alpha\log\beta - n\log\Gamma(\alpha) + (\alpha - 1)\, n\,\overline{\log x} - \beta\, n\bar{x}\]

With numbers:

\[\ell(\alpha, \beta) = 1000\left[\alpha\log\beta - \log\Gamma(\alpha) + (\alpha - 1)(5.078763) - 160.7419\,\beta\right]\]

loglik_gamma <- function(a, b) n * (a * log(b) - lgamma(a) + (a - 1) * mlog - b * xbar)
c(formula = loglik_gamma(a_mm, b_mm), dgamma_sum = sum(dgamma(x, shape = a_mm, rate = b_mm, log = TRUE)))
##    formula dgamma_sum 
##  -3408.705  -3408.705

Step 2. Solve for the rate

\[\frac{\partial \ell}{\partial \beta} = \frac{n\alpha}{\beta} - n\bar{x} = 0 \;\Longrightarrow\; \frac{\alpha}{\beta} = \bar{x} \;\Longrightarrow\; \hat{\beta}_{\text{MLE}} = \frac{\hat{\alpha}_{\text{MLE}}}{\bar{x}}\]

So the maximum likelihood estimation (MLE) fit always matches the sample mean exactly: \(\hat\alpha/\hat\beta = \bar{x}\).

Step 3. The shape equation has no formula

With \(\psi(\alpha) = \frac{d}{d\alpha}\log\Gamma(\alpha)\):

\[\frac{\partial \ell}{\partial \alpha} = n\log\beta - n\,\psi(\alpha) + n\,\overline{\log x} = 0\]

Substitute \(\beta = \alpha/\bar{x}\) from Step 2:

\[\log\alpha - \log\bar{x} - \psi(\alpha) + \overline{\log x} = 0 \;\Longrightarrow\; \log\alpha - \psi(\alpha) = \log\bar{x} - \overline{\log x}\]

With numbers, the right side is \(5.079800 - 5.078763 = 0.001037\):

\[\log\alpha - \psi(\alpha) = 0.001037\]

When we take the slope with respect to α, α ends up both inside a log and inside another function built from the Gamma function. There is no way to rearrange the equation to get α by itself, so the computer has to solve it numerically (textbook Section 8.4.3). R calls that function digamma(). All you need to know is that it has no simple formula. One clue helps: the right side of the equation is the gap from Part 0, 0.001037. A gap that small means α must be large, around 482.

gap <- log(xbar) - mlog
c(rhs = gap, lhs_at_a_mm = log(a_mm) - digamma(a_mm), rough_alpha = 1 / (2 * gap))
##          rhs  lhs_at_a_mm  rough_alpha 
## 1.036935e-03 1.037297e-03 4.821903e+02

Step 4. Let the computer search

Negative log-likelihood on the log scale, \(p = (\log\alpha, \log\beta)\):

\[\text{nll}(p) = -\sum_{i=1}^{n}\log f\!\left(x_i;\ e^{p_1}, e^{p_2}\right)\]

Start: \(p_0 = (\log 482.19, \log 2.9998)\).

R’s optim() finds the best α and β by trial and improvement. It starts at the method of moments (MM) values, tries nearby values, keeps whichever makes the data more likely, and stops when the improvements get very small. With R’s default stopping setting it quit too early, at α = 482.47. Near the top, the likelihood barely changes along the line β = α/x̄, so the improvements looked small before it had reached the peak. Telling it to keep going longer (reltol = 1e-12 in the code) gave α = 482.36 and β = 3.0008. The search works on log α and log β so it can never try a negative value. Solving the Step 3 equation a second way, with uniroot(), gives the same α = 482.36.

\[\hat{\alpha}_{\text{MLE}} = 482.36, \qquad \hat{\beta}_{\text{MLE}} = 3.0008, \qquad \texttt{fit\$convergence} = 0\]

The maximized log-likelihood is \(-3408.7050\), slightly higher than the Normal’s \(-3409.0706\).

fit_gamma_mle <- function(v, start, reltol = 1e-12) {
  nll <- function(p) -sum(dgamma(v, shape = exp(p[1]), rate = exp(p[2]), log = TRUE))
  optim(log(start), nll, control = list(reltol = reltol))
}
fit_default <- fit_gamma_mle(x, c(a_mm, b_mm), reltol = 1e-8)
fit <- fit_gamma_mle(x, c(a_mm, b_mm))
a_mle <- exp(fit$par[1]); b_mle <- exp(fit$par[2])
a_exact <- uniroot(function(a) log(a) - digamma(a) - gap, c(1, 1e4), tol = 1e-12)$root
c(default_a = exp(fit_default$par[1]), default_b = exp(fit_default$par[2]),
  a_mle = a_mle, b_mle = b_mle, convergence = fit$convergence,
  a_uniroot = a_exact, loglik = -fit$value)
##    default_a    default_b        a_mle        b_mle  convergence    a_uniroot 
##   482.465629     3.001518   482.356598     3.000814     0.000000   482.356894 
##       loglik 
## -3408.704954

Step 5. Confirm the answer

\[\frac{\hat\alpha}{\bar{x}} = \frac{482.36}{160.7419} = 3.000814 \qquad \text{vs} \qquad \hat\beta = 3.000814\]

They agree to 7 significant digits, so the optimizer did land on the ridge \(\beta = \alpha/\bar{x}\) that Step 2 requires.

c(b_mle = b_mle, a_mle_over_xbar = a_mle / xbar, model_mean = a_mle / b_mle)
##           b_mle a_mle_over_xbar      model_mean 
##        3.000814        3.000814      160.741896

How the two Gamma fits compare

shape \(\hat\alpha\) rate \(\hat\beta\) mean \(\hat\alpha/\hat\beta\)
method of moments (MM) 482.19 2.9998 160.7419
maximum likelihood estimation (MLE) 482.36 3.0008 160.7419
change 0.17 (0.035 percent) 0.0010 (0.035 percent) 0.0000

The shape moved more in absolute terms (up 0.17 vs up 0.0010); in relative terms both rose by the same 0.035 percent, because the rate is tied to the shape through \(\beta = \alpha/\bar{x}\). The mean does not change: both methods force \(\hat\alpha/\hat\beta = \bar{x}\). Only the implied variance \(\hat\alpha/\hat\beta^2 = \bar{x}/\hat\beta\) changes, from 53.5848 to 53.5661. The maximum likelihood estimation (MLE) variance sits between the \(n - 1\) variance 53.5848 and the \(1/n\) variance 53.5312.

Part 4: Exponential, by method of moments (MM) and maximum likelihood estimation (MLE)

The Exponential is a one-dial model built for waiting times. Using it for height is a stress test: the homework includes it on purpose so you can see what a wrong family does.

Step 1. Method of moments (MM)

Model mean \(= 1/\lambda\). Set equal to the sample mean:

\[\frac{1}{\lambda} = \bar{x} \;\Longrightarrow\; \hat\lambda_{\text{MM}} = \frac{1}{\bar{x}} = \frac{1}{160.7419} = 0.0062212 \text{ per cm}\]

lam_mm <- 1 / xbar
lam_mm
## [1] 0.006221153

Step 2. Maximum likelihood estimation (MLE)

Density \(f(x; \lambda) = \lambda e^{-\lambda x}\). Log of one term: \(\log\lambda - \lambda x_i\). Sum over \(n\):

\[\ell(\lambda) = n\log\lambda - \lambda\sum_{i=1}^{n} x_i\]

\[\frac{d\ell}{d\lambda} = \frac{n}{\lambda} - \sum_{i=1}^{n} x_i = 0 \;\Longrightarrow\; \hat\lambda = \frac{n}{\sum x_i} = \frac{1}{\bar{x}}\]

\[\frac{d^2\ell}{d\lambda^2} = -\frac{n}{\lambda^2} < 0 \;\text{(maximum)}\]

\[\hat\lambda_{\text{MLE}} = \frac{1000}{160741.9} = \frac{1}{160.7419} = 0.0062212 \text{ per cm}\]

Same as the method of moments (MM). Maximized log-likelihood: \(\ell(\hat\lambda) = -6079.8000\), far below the Normal’s \(-3409.0706\).

lam_mle <- n / sum(x)
num <- optimize(function(l) n * log(l) - l * sum(x), c(1e-5, 1), maximum = TRUE, tol = 1e-12)
c(lam_mle = lam_mle, optimize_check = num$maximum, loglik = n * log(lam_mle) - lam_mle * sum(x))
##        lam_mle optimize_check         loglik 
##   6.221153e-03   6.221153e-03  -6.079800e+03

Step 3. Compare the fit to the data

value (cm)
fitted mean \(1/\hat\lambda\) 160.7419
fitted peak (mode) 0
sample mean \(\bar{x}\) 160.7419
fitted standard deviation \(1/\hat\lambda\) 160.7419
sample standard deviation \(s\) 7.3202

The Exponential matches the mean by construction, but it puts its peak at 0 cm and its standard deviation equal to its mean, about 22 times the real spread.

What to expect from the plot

The histogram is a tall narrow hump from about 139 to 191 cm reaching a density near 0.055. The Exponential curve starts at height \(\hat\lambda = 0.0062212\) at 0 cm and decays slowly, reaching about \(0.00229\) at 160.7 cm. So the curve is a nearly flat, low line that is highest where there are no women and far below the bars where the women are.

Part 5: The fit visuals

A fitted distribution is only useful if it looks like the data. We check all six fits three ways, because each plot catches a different kind of mismatch.

Six fits: Normal, Gamma, Exponential, each by method of moments (MM) and by maximum likelihood estimation (MLE).

Step 1. Density over histogram

Six panels below, one per fit. Normal and Gamma panels use 135 to 195 cm; Exponential panels use 0 to 250 cm so the curve’s peak at 0 is visible. Histogram: 30 bins on the density scale (freq = FALSE).

fits <- list(
  list(fam = "Normal", meth = "method of moments (MM)",
       d = function(v) dnorm(v, mu_mm, sigma_mm), p = function(v) pnorm(v, mu_mm, sigma_mm),
       q = function(pr) qnorm(pr, mu_mm, sigma_mm), lim = c(135, 195)),
  list(fam = "Normal", meth = "maximum likelihood estimation (MLE)",
       d = function(v) dnorm(v, mu_mle, sigma_mle), p = function(v) pnorm(v, mu_mle, sigma_mle),
       q = function(pr) qnorm(pr, mu_mle, sigma_mle), lim = c(135, 195)),
  list(fam = "Gamma", meth = "method of moments (MM)",
       d = function(v) dgamma(v, shape = a_mm, rate = b_mm), p = function(v) pgamma(v, shape = a_mm, rate = b_mm),
       q = function(pr) qgamma(pr, shape = a_mm, rate = b_mm), lim = c(135, 195)),
  list(fam = "Gamma", meth = "maximum likelihood estimation (MLE)",
       d = function(v) dgamma(v, shape = a_mle, rate = b_mle), p = function(v) pgamma(v, shape = a_mle, rate = b_mle),
       q = function(pr) qgamma(pr, shape = a_mle, rate = b_mle), lim = c(135, 195)),
  list(fam = "Exponential", meth = "method of moments (MM)",
       d = function(v) dexp(v, rate = lam_mm), p = function(v) pexp(v, rate = lam_mm),
       q = function(pr) qexp(pr, rate = lam_mm), lim = c(0, 250)),
  list(fam = "Exponential", meth = "maximum likelihood estimation (MLE)",
       d = function(v) dexp(v, rate = lam_mle), p = function(v) pexp(v, rate = lam_mle),
       q = function(pr) qexp(pr, rate = lam_mle), lim = c(0, 250)))
ttl <- function(f) paste0(f$fam, "\n", f$meth)

dark(); par(mfrow = c(2, 3))
for (f in fits) {
  hist(x, freq = FALSE, breaks = 30, xlim = f$lim, ylim = c(0, 0.065), col = "gray30", border = "gray60",
       main = ttl(f), xlab = "height (cm)")
  curve(f$d(x), from = f$lim[1], to = f$lim[2], n = 1001, add = TRUE, col = "deepskyblue", lwd = 3)
}

Step 2. Cumulative curves

The Kolmogorov-Smirnov (KS) distance is the largest vertical gap between the fitted curve and the data’s curve. 0 means perfect agreement, and 1 is the worst possible.

dark(); par(mfrow = c(2, 3))
Fn <- ecdf(x); xs <- sort(x)
ks_gap <- function(f) max(abs(c(Fn(xs) - f$p(xs), Fn(xs) - 1 / n - f$p(xs))))
for (f in fits) {
  plot(Fn, do.points = FALSE, verticals = TRUE, xlim = f$lim, col = ink, lwd = 2,
       main = ttl(f), xlab = "height (cm)", ylab = "cumulative probability")
  curve(f$p(x), from = f$lim[1], to = f$lim[2], n = 1001, add = TRUE, col = "deepskyblue", lwd = 3, lty = 2)
}

sapply(fits, ks_gap)
## [1] 0.01545283 0.01533868 0.01689911 0.01686009 0.58170101 0.58170101

Step 3. Quantile-quantile (QQ) plot

A quantile-quantile (QQ) plot pairs each sorted height with where the fitted distribution says that height should fall. Points on the 45 degree line mean the two agree.

dark(); par(mfrow = c(2, 3))
for (f in fits) {
  model_q <- f$q(ppoints(n))
  plot(model_q, sort(x), pch = 16, cex = 0.6, col = "deepskyblue", main = ttl(f),
       xlab = "model quantile (cm)", ylab = "sample quantile (cm)",
       xlim = range(model_q, x), ylim = range(model_q, x))
  abline(0, 1, col = "orange", lwd = 2)
}

range(fits[[5]]$q(ppoints(n)))
## [1] 8.039105e-02 1.221784e+03

Step 4. All six fits at a glance

Fit Probability density function (PDF) over histogram Cumulative distribution function (CDF) over empirical cumulative distribution function (ECDF) Quantile-quantile (QQ) plot
Normal, method of moments (MM) Center and spread match; tails follow the bars. Hugs the steps everywhere; largest gap 0.0155. On the line from 145 to 175 cm; tallest women above the line (upper tail slightly longer than the model).
Normal, maximum likelihood estimation (MLE) Indistinguishable from the method of moments (MM) fit. Hugs the steps; largest gap 0.0153. Indistinguishable from the method of moments (MM) fit.
Gamma, method of moments (MM) Indistinguishable from the Normal. Hugs the steps; largest gap 0.0169. On the line; extremes slightly closer to it than the Normal’s, a few points near 142 cm just below.
Gamma, maximum likelihood estimation (MLE) Indistinguishable from the Gamma method of moments (MM) fit. Hugs the steps; largest gap 0.0169. Same as the Gamma method of moments (MM) fit.
Exponential, method of moments (MM) Wrong shape: peak at 0 cm, curve far below the bars. Wrong everywhere; largest gap 0.5817. Strong curve far off the line; center, spread, and tails all wrong.
Exponential, maximum likelihood estimation (MLE) Identical to the method of moments (MM) fit. Identical; largest gap 0.5817. Identical.

Part 6: The median and its sampling distribution

The median is the height that splits the women in half. Each fitted distribution has its own median. Because each fit came from one sample, a different sample would give a slightly different median. Here we compute each median and measure how much it could vary.

Step 1. Median of each fitted distribution

The model median \(m\) solves \(F(m; \hat\theta) = 0.5\), which is what the q-function returns at 0.5.

Fit Model median (cm) median(x) (cm) Model minus data
Normal, method of moments (MM) 160.7419 160.8 -0.0581
Normal, maximum likelihood estimation (MLE) 160.7419 160.8 -0.0581
Gamma, method of moments (MM) 160.6308 160.8 -0.1692
Gamma, maximum likelihood estimation (MLE) 160.6308 160.8 -0.1692
Exponential, method of moments (MM) 111.4178 160.8 -49.3822
Exponential, maximum likelihood estimation (MLE) 111.4178 160.8 -49.3822 
model_medians <- c(normal_mm = qnorm(0.5, mu_mm, sigma_mm),
                   normal_mle = qnorm(0.5, mu_mle, sigma_mle),
                   gamma_mm = qgamma(0.5, shape = a_mm, rate = b_mm),
                   gamma_mle = qgamma(0.5, shape = a_mle, rate = b_mle),
                   exp_mm = qexp(0.5, rate = lam_mm),
                   exp_mle = qexp(0.5, rate = lam_mle))
round(cbind(model = model_medians, data = median(x), model_minus_data = model_medians - median(x)), 4)
##               model  data model_minus_data
## normal_mm  160.7419 160.8          -0.0581
## normal_mle 160.7419 160.8          -0.0581
## gamma_mm   160.6308 160.8          -0.1692
## gamma_mle  160.6308 160.8          -0.1692
## exp_mm     111.4178 160.8         -49.3822
## exp_mle    111.4178 160.8         -49.3822

Step 2. Which medians land near the data

Normal and Gamma medians land within 0.2 cm of median(x) \(= 160.8\) (Normal \(-0.0581\), Gamma \(-0.1692\) and \(-0.1692\)). The Gamma median sits a hair left of its mean because of its slight right skew. These are the fits whose curves traced the histogram and whose quantile-quantile (QQ) points sat on the line. The Exponential median is \(111.4178\) cm, \(49.38\) cm below the data median, matching the quantile-quantile (QQ) plot where every model percentile was in the wrong place.

Step 3. Simulate one fit

How much could the median change with a different sample? To find out, we simulate. We treat the fitted distribution as if it were the truth, draw a fake sample of 1,000 heights from it, refit it the same way, and record its median. We repeat that 5,000 times. The spread of those 5,000 medians is the median’s sampling distribution (textbook Section 7.2), and its standard deviation is the standard error.

Fit: Normal, method of moments (MM). \(B = 5000\), seed set once at the top (set.seed(6110)), no other random calls before this point.

For \(b = 1, \dots, B\):

  1. Draw \(x^*_1, \dots, x^*_{1000} \sim \text{Normal}(160.7419, 7.3202)\).
  2. Re-estimate: \(\hat\mu^* = \bar{x}^*\), \(\hat\sigma^* = \sqrt{\frac{1}{n-1}\sum (x^*_i - \bar{x}^*)^2}\) (sd()).
  3. Model median: \(m^*_b\) = qnorm(0.5, mu_star, sigma_star), which equals \(\hat\mu^*\).

Result: standard error \(= \operatorname{sd}(m^*) = 0.2334\) cm; middle 95 percent from 160.2878 to 161.2153 cm.

Sanity check: the model median of a Normal fit equals \(\bar{x}^*\), whose standard error is \(\hat\sigma/\sqrt{n} = 7.3202/\sqrt{1000} = 0.2315\) cm. The simulation gives 0.2334.

B <- 5000
refit_normal_mm <- function(s) qnorm(0.5, mean(s), sd(s))                        # divide by n - 1
refit_normal_mle <- function(s) { m <- mean(s); qnorm(0.5, m, sqrt(mean((s - m)^2))) }  # divide by n
boot_normal_mm <- replicate(B, refit_normal_mm(rnorm(n, mu_mm, sigma_mm)))
c(se = sd(boot_normal_mm), quantile(boot_normal_mm, c(0.025, 0.975)), theory_se = sigma_mm / sqrt(n))
##          se        2.5%       97.5%   theory_se 
##   0.2334025 160.2877582 161.2152591   0.2314838

Step 4. Simulate all six

Each fit is refit with its own method inside the loop: Normal maximum likelihood estimation (MLE) with the \(1/n\) standard deviation; Gamma method of moments (MM) with \(\bar{x}^{*2}/s^{*2}\) and \(\bar{x}^*/s^{*2}\) using var(); Gamma maximum likelihood estimation (MLE) with optim() started at that fake sample’s method of moments (MM) values (all \(B\) runs reported convergence code 0: 0 failures); Exponential with \(1/\bar{x}^*\).

Fit Model median Simulated standard error 2.5th percentile 97.5th percentile Contains median(x) = 160.8?
Normal, method of moments (MM) 160.7419 0.2334 160.2878 161.2153 yes
Normal, maximum likelihood estimation (MLE) 160.7419 0.2303 160.2941 161.1987 yes
Gamma, method of moments (MM) 160.6308 0.2338 160.1655 161.0856 yes
Gamma, maximum likelihood estimation (MLE) 160.6308 0.2268 160.1913 161.0746 yes
Exponential, method of moments (MM) 111.4178 3.5122 104.6399 118.5099 no
Exponential, maximum likelihood estimation (MLE) 111.4178 3.4731 104.7166 118.5363 no

The Normal and Gamma histograms are narrow bells centered on their model medians with standard error about 0.23 cm. The Exponential histogram is centered near 111 cm with standard error about 3.5122 cm, about 15 times wider, because the Exponential’s spread is tied to its mean.

refit_gamma_mm <- function(s) { m <- mean(s); vs <- var(s); qgamma(0.5, shape = m^2 / vs, rate = m / vs) }
nonconv <- 0
refit_gamma_mle <- function(s) {
  m <- mean(s); vs <- var(s)
  f <- fit_gamma_mle(s, c(m^2 / vs, m / vs))   # start at this sample's method of moments (MM) values
  if (f$convergence != 0) nonconv <<- nonconv + 1
  qgamma(0.5, shape = exp(f$par[1]), rate = exp(f$par[2]))
}
refit_exp <- function(s) qexp(0.5, rate = 1 / mean(s))

boot_normal_mle <- replicate(B, refit_normal_mle(rnorm(n, mu_mle, sigma_mle)))
boot_gamma_mm   <- replicate(B, refit_gamma_mm(rgamma(n, shape = a_mm, rate = b_mm)))
boot_gamma_mle  <- replicate(B, refit_gamma_mle(rgamma(n, shape = a_mle, rate = b_mle)))
boot_exp_mm     <- replicate(B, refit_exp(rexp(n, rate = lam_mm)))
boot_exp_mle    <- replicate(B, refit_exp(rexp(n, rate = lam_mle)))

boots <- list(normal_mm = boot_normal_mm, normal_mle = boot_normal_mle, gamma_mm = boot_gamma_mm,
              gamma_mle = boot_gamma_mle, exp_mm = boot_exp_mm, exp_mle = boot_exp_mle)
boot_table <- t(sapply(boots, function(v) c(se = sd(v), quantile(v, c(0.025, 0.975)))))
print(round(cbind(model_median = model_medians, boot_table), 4))
##            model_median     se     2.5%    97.5%
## normal_mm      160.7419 0.2334 160.2878 161.2153
## normal_mle     160.7419 0.2303 160.2941 161.1987
## gamma_mm       160.6308 0.2338 160.1655 161.0856
## gamma_mle      160.6308 0.2268 160.1913 161.0746
## exp_mm         111.4178 3.5122 104.6399 118.5099
## exp_mle        111.4178 3.4731 104.7166 118.5363
cat("Gamma maximum likelihood refits that did not converge:", nonconv, "\n")
## Gamma maximum likelihood refits that did not converge: 0
dark(); par(mfrow = c(2, 3))
for (i in seq_along(boots)) {
  hist(boots[[i]], breaks = 40, freq = FALSE, col = "gray30", border = "gray60",
       main = ttl(fits[[i]]), xlab = "simulated model median (cm)")
  abline(v = model_medians[i], col = "deepskyblue", lwd = 3)
  abline(v = quantile(boots[[i]], c(0.025, 0.975)), col = "orange", lwd = 2, lty = 2)
  abline(v = median(x), col = "red", lwd = 2, lty = 3)
}

Does a narrow interval mean the model is right?

No. The Exponential’s 95 percent simulated interval, 104.6399 to 118.5099 cm, is only about 13.9 cm wide, yet even its upper end sits about 42 cm below the data median 160.8 cm and does not contain it. The simulation treats the fitted model as the truth, so it only measures how much the estimate would wobble if the model were right. It says nothing about whether the family is right. Model checking (Part 5) answers that question.

Checklist

Every item on the assignment’s checklist, for both methods and all three distributions, and where to find it in this tutorial.

Item Normal, MM Normal, MLE Gamma, MM Gamma, MLE Exponential, MM Exponential, MLE
Estimates of parameters μ = 160.7419, σ = 7.3202 μ = 160.7419, σ = 7.3165 α = 482.19, β = 2.9998 α = 482.36, β = 3.0008 λ = 0.006221 λ = 0.006221
Estimated PDF over histogram ✓ Part 5, Step 1 ✓ Part 5, Step 1 ✓ Part 5, Step 1 ✓ Part 5, Step 1 ✓ Part 5, Step 1 ✓ Part 5, Step 1
Estimated CDF over ECDF ✓ Part 5, Step 2 ✓ Part 5, Step 2 ✓ Part 5, Step 2 ✓ Part 5, Step 2 ✓ Part 5, Step 2 ✓ Part 5, Step 2
QQ plot (sample vs. estimated) ✓ Part 5, Step 3 ✓ Part 5, Step 3 ✓ Part 5, Step 3 ✓ Part 5, Step 3 ✓ Part 5, Step 3 ✓ Part 5, Step 3
Estimated median (simulated standard error) 160.7419 (0.2334) 160.7419 (0.2303) 160.6308 (0.2338) 160.6308 (0.2268) 111.4178 (3.5122) 111.4178 (3.4731)

MM = method of moments; MLE = maximum likelihood estimation. The real median of the 1,000 heights is 160.8 cm.

Part 7: Compare, contrast, and take-home messages

Step back and compare. Which mattered more: the method or the choice of distribution?

Step 1. Method versus method

Step 2. Distribution versus distribution

  1. Gamma and Normal, essentially tied. Gamma has the higher log-likelihood (-3408.7050 vs -3409.0706, a difference of 0.3657); Normal has the smaller largest cumulative distribution function (CDF) gap (0.0153 vs 0.0169) and its median is closer to median(x) (-0.0581 vs -0.1692). With a Gamma shape near 480 the two curves are visually identical. For a newcomer, the Normal is the simpler choice, because its two parameters are the mean and the standard deviation directly.
  2. Exponential, a distant last. Log-likelihood -6079.8000, largest cumulative distribution function (CDF) gap 0.5817, model median off by 49.38 cm.

Step 3. Which choice mattered more

The family. Changing the method moved the Gamma model median by 0.00003 cm. Changing the family from Normal to Exponential moved the model median by 49.32 cm, from 160.7419 to 111.4178. The family change was over a million times larger.

Step 4. Assumptions

  1. The 1,000 heights are independent draws from one common distribution. The National Health and Nutrition Examination Survey (NHANES) uses a complex survey design with sampling weights; these are ignored here.
  2. The sample is the first 1,000 rows after filtering to women aged 18 or older, as the homework specifies. This is not a random subsample, so row order could matter.
  3. Units are centimeters. No missing values were present, so nothing was dropped.
  4. Method of moments (MM) variance divides by \(n - 1\) (var() \(= 53.5848\)). The \(1/n\) version (\(53.5312\)) appears only as the Normal maximum likelihood estimation (MLE).
  5. Simulation: 5,000 fake samples of 1,000 heights for each fit, with set.seed(6110) set once at the top, run in the order Normal method of moments (MM), Normal maximum likelihood estimation (MLE), Gamma method of moments (MM), Gamma maximum likelihood estimation (MLE), Exponential method of moments (MM), Exponential maximum likelihood estimation (MLE).
  6. Gamma maximum likelihood estimation (MLE): found by optim() searching over log α and log β, starting at the method of moments (MM) values, with a tighter stopping setting than R’s default because the default quit early. Confirmed by solving the equation a second way with uniroot().
  7. The Gamma is parameterized by shape and rate; every R call names shape = and rate =.

Step 5. Take-home messages

  1. Pick the family first, then worry about the method. For these heights, switching between method of moments (MM) and maximum likelihood estimation (MLE) changed the Gamma median by 0.00003 cm. Switching from a Normal to an Exponential changed it by 49.32 cm.
  2. Method of moments (MM) is algebra, maximum likelihood estimation (MLE) is optimization. Method of moments (MM) sets model averages equal to data averages and solves. Maximum likelihood estimation (MLE) finds the parameter values that make the observed data most probable; sometimes calculus gives a formula (Normal, Exponential), sometimes a computer has to search (Gamma shape). For the Exponential the two methods give identical answers; for the Normal they differ only through the variance divisor.
  3. Always plot the fit three ways. A density curve over a histogram, a cumulative distribution function (CDF) over the empirical cumulative distribution function (ECDF), and a quantile-quantile (QQ) plot each catch different problems. The Exponential matched the sample mean exactly and still failed all three.
  4. A confident answer can still be wrong. The Exponential’s simulated interval for the median was narrow and was centered about 49 cm below the real median. The simulation measures wobble under the model, not whether the model is right.
  5. Check the software defaults. R’s Gamma functions use rate unless you say otherwise, var() divides by n − 1, and optim() can stop searching before it reaches the peak. Name every argument, state your divisor, and confirm the computer’s answer a second way.

↑ Back to top