5. Selecting Input Probability Distributions

How to choose the input distributions a simulation feeds on: assess sample independence, hypothesize a distribution family, estimate its parameters by maximum likelihood, and judge the fit with chi-square and K-S tests.

flowchart LR
    RNG[RNG] -->|"U(0,1)"| T[transformations]
    F(["f / F (input distribution)"]) --> T
    T -->|RVs| P((simulation program))

Chapter 3 built the RNG and Chapter 4 the transformations. This chapter is about the remaining box: which input distribution $f$ / $F$ to use.

E.g. for the $M/M/1$ queue we simply assumed exponential interarrival and service times. Q:

1. How to select input distributions?
2. Is a specified input distribution a good approximation of the actual process?

flowchart LR
    D[Data] --> I[input distributions]

In practice, an input distribution is specified based on data collected from the actual process by observing an input r.v. of interest.

Three approaches

1. Use the actual data directly (e.g. customer service times): store the data in a “lookup” table and read the values in order during the simulation (trace-driven simulation).
   1. Disadvantage: this approach is very expensive and time-consuming (a large amount of data must be stored and read).
   2. Disadvantage: it can only reproduce what has happened historically.
   3. Use case: typically used in failure analysis, and to validate a model of an existing system (feed the model the same historical inputs and compare its output with the real system’s).
2. Use the data to fit an empirical distribution (frequency distribution, histogram), and then sample from the empirical distribution during the simulation.
   1. Disadvantage: it is not possible to generate values outside the range of the observed data $X_1, \dots, X_n$, i.e. outside $\left[\min_i X_i,\ \max_i X_i\right]$.
3. Use the data to fit a theoretical distribution (a parameterized distribution such as $\text{Exp}(\lambda)$, $N(\mu, \sigma^2)$, $U(a, b)$); use statistical methods to determine the quality of fit.
   1. Advantage: compact representation of the data (a few parameters instead of $n$ numbers).
   2. Advantage: flexible → “what-if” questions (e.g. what if service becomes 10% faster? Just change a parameter).
   3. Advantage: it can generate values outside the observed range, and it smooths out the irregularities of a small sample.

Use the third approach whenever possible, i.e. whenever some theoretical distribution fits the data well.

Assumption: to fit a distribution, the data should be i.i.d. observations, i.e. $X_1, X_2, \dots, X_n$ are i.i.d. from an underlying distribution. (§5.1 checks the independence part.)

5.1. Techniques for assessing sample independence

Given $X_1, X_2, \dots, X_n$, in the order they were observed.

• Correlation plot

The correlation between two observations $j$ apart (lag $j$) is

$$\rho_j = \operatorname{corr}(X_i, X_{i+j}) = \frac{\operatorname{cov}(X_i, X_{i+j})}{\sqrt{\operatorname{Var}(X_i)}\,\sqrt{\operatorname{Var}(X_{i+j})}} .$$

If $\rho_j \ne 0$ for some $j$, then the $X_i$’s are not independent.

From the data, estimate $\rho_j$ by

$$\hat\rho_j = \frac{\hat C_j}{S^2(n)}, \qquad \hat C_j = \frac{1}{n-j}\sum_{i=1}^{n-j}\big(X_i - \bar X(n)\big)\big(X_{i+j} - \bar X(n)\big),$$

and plot $\hat\rho_j$ against $j = 1, 2, \dots$: this is the correlation plot. Even for independent data $\hat\rho_j$ is never exactly 0, but it should stay close to 0 (a common rule of thumb: within $\pm 2/\sqrt n$). Values clearly away from 0 indicate dependence. Since a correlation of 0 does not prove independence (§2.3), the plot can only detect dependence, not prove independence.

• Scatter diagram

Group the observations into pairs $(X_i, X_{i+1})$, $i = 1, 2, \dots, n-1$, and plot them as points.

If $X_i$ and $X_{i+1}$ are positively correlated, then the points will approximately lie along a line with a positive slope (negatively correlated: negative slope). If they are independent, the points are scattered with no pattern.

• Runs test

The runs-up test of Chapter 3 (§3.4.3). It uses only the order of the observations (which one is larger than which), so it works for any continuous data, not just $U(0,1)$ numbers.

Three steps to fit a theoretical distribution

1. Hypothesize a family of distributions (§5.4).
2. Estimate the parameters of the distribution (maximum likelihood estimator, MLE).
3. Determine the quality of fit (validation).

5.2. Review of distributions and their parameters

The types of parameters:

• Location parameter $\gamma$

It determines the location of the distribution (where it sits on the $x$-axis). Changing $\gamma$ shifts the distribution but does not change its shape.

If $X$ has location parameter $0$, then $Y = X + \gamma$ has location parameter $\gamma$. E.g. for $N(\mu, \sigma^2)$, $\mu$ is the location parameter: if $X \sim N(0, \sigma^2)$, then $X + \mu \sim N(\mu, \sigma^2)$.

• Scale parameter $\beta$

It determines the scale / statistical dispersion of the distribution. Changing $\beta$ expands / compresses the distribution, but the distribution maintains its basic shape. If $X$ has scale parameter $1$, then $Y = \beta X$ has scale parameter $\beta$.

E.g.

$N(\mu, \sigma^2)$	$\text{Exp}(\lambda)$	$U(a, b)$
$\sigma$ is the scale parameter: if $X \sim N(0,1)$, then $\sigma X \sim N(0, \sigma^2)$. ($\mu$ is the location parameter.)	$\operatorname{Var}(X) = 1/\lambda^2$. The scale parameter is $\beta = 1/\lambda$, the mean: if $X \sim \text{Exp}(1)$, then $\beta X$ is exponential with mean $\beta$. $\lambda$ itself is a rate, so a large $\lambda$ compresses the distribution.	$f(x) = \dfrac{1}{b-a}$. Let $L = b - a$ and reparameterize $U(a, b)$ by $(a, L)$: $a$ is the location and $L$ the scale parameter, since $a + L\,U \sim U(a, a+L)$ for $U \sim U(0,1)$. The density height is $1/L$.

• Shape parameter $\alpha$

It determines the basic form or shape of the distribution, not just its position and spread.

Ex. Weibull distribution ($\alpha$: shape, $\beta$: scale):

$$f(x) = \alpha\beta^{-\alpha}x^{\alpha-1}e^{-(x/\beta)^{\alpha}}, \qquad x > 0 .$$

$\alpha < 1$: decreasing density, with $f(x) \to \infty$ as $x \to 0$. $\alpha = 1$: the exponential distribution with mean $\beta$. $\alpha > 1$: a single hump, with $f(0) = 0$. (The normal and exponential families have no shape parameter: every normal density has the same bell shape.)

5.3. Fit an empirical distribution

• Used when there is enough data and theoretical distributions do not work well.

Given $n$ observations $X_1, \dots, X_n$, sort them: $X_{(1)} < X_{(2)} < \dots < X_{(n)}$ (order statistics). $F(X_{(i)}) = i/n$.

Empirical c.d.f.:

$$F(x) = \frac{\text{number of } X_i\text{'s} \le x}{n} .$$

For continuous data, use the interpolated (piecewise linear) version:

$$\tilde F(x) = \begin{cases} 0, & x < X_{(1)}, \\[4pt] \dfrac{i-1}{n-1} + \dfrac{x - X_{(i)}}{(n-1)\big(X_{(i+1)} - X_{(i)}\big)}, & X_{(i)} \le x < X_{(i+1)}, \quad i = 1, \dots, n-1, \\[10pt] 1, & x \ge X_{(n)} . \end{cases}$$

$$\tilde F(X_{(i)}) = \frac{i-1}{n-1}, \qquad F(X_{(i)}) = \frac{\text{number of } X_k\text{'s} \le X_{(i)}}{n} = \frac{i}{n} .$$

Drawbacks

1. The mean of $\tilde F$ $\ne \bar X$. $\tilde F$ spreads probability $\frac{1}{n-1}$ uniformly over each interval $[X_{(i)}, X_{(i+1)}]$, so its mean is $\frac{1}{n-1}\sum_{i=1}^{n-1}\frac{X_{(i)} + X_{(i+1)}}{2}$, which in general differs from $\bar X$. E.g. for the data $1.0, 2.2, 2.8, 4.5, 5.6$: mean of $\tilde F$ $= 3.20$, $\bar X = 3.22$.
2. Random values generated from $\tilde F$ can never be smaller than $X_{(1)}$ or larger than $X_{(n)}$.

Sampling from $\tilde F$ (inverse transform, §4.2.1): generate $u \sim U(0,1)$, let $P = (n-1)u$ and $I = \lfloor P \rfloor + 1$, and return $X = X_{(I)} + (P - I + 1)\big(X_{(I+1)} - X_{(I)}\big)$.

5.4. Step I: Hypothesize a distribution family

• Prior knowledge of the random variable of interest
  • Service times should never be fitted using a normal distribution: a normal r.v. can be negative, a service time cannot (and service times are usually right-skewed, while the normal is symmetric).
  • Customers arrive at a bank one at a time (independently, at a roughly constant rate) ⇒ the arrival process is a Poisson process, so the interarrival times are exponential (§4.6).
• Other heuristics
  • Summary statistics: mean, variance, min, max, mode. Also the median (mean ≈ median suggests a symmetric distribution) and the coefficient of variation $\sigma/\mu$ (equal to 1 for the exponential).
  • Histogram (choose an appropriate number of subintervals: too few hides the shape, too many makes it ragged; a common starting point is $k \approx 1 + \log_2 n$, then try a few others). The shape of the histogram estimates the shape of the density.
  • Quantile summaries and boxplot.

Given $q \in (0, 1)$, the $q$-quantile is the value $x_q$ such that $P(X \le x_q) \ge q$ and $P(X \ge x_q) \ge 1 - q$. (For continuous $F$: $F(x_q) = q$, i.e. $x_q = F^{-1}(q)$. $x_{0.5}$ is the median; $x_{0.25}$ and $x_{0.75}$ are the quartiles.)

For a given sample $X_1, X_2, \dots, X_n$ with $X_{(1)} \le X_{(2)} \le \dots \le X_{(n)}$, the sample $100q\%$ quantile is

$$\hat x_q = \begin{cases} X_{(\lceil nq \rceil)}, & \text{if } nq \text{ is not an integer}, \\[6pt] \dfrac{X_{(nq)} + X_{(nq+1)}}{2}, & \text{if } nq \text{ is an integer}. \end{cases}$$

Ex. $n = 10$ service times, sorted: $0.2,\ 0.5,\ 0.7,\ 1.1,\ 1.4,\ 1.8,\ 2.3,\ 3.0,\ 4.2,\ 6.5$.

• $q = 0.25$: $nq = 2.5$ is not an integer ⇒ $\hat x_{0.25} = X_{(3)} = 0.7$.
• $q = 0.5$: $nq = 5$ is an integer ⇒ $\hat x_{0.5} = \dfrac{X_{(5)} + X_{(6)}}{2} = \dfrac{1.4 + 1.8}{2} = 1.6$.
• $q = 0.75$: $nq = 7.5$ ⇒ $\hat x_{0.75} = X_{(8)} = 3.0$.

Boxplot

Draws the five numbers $x_{\min},\ \hat x_{0.25},\ \hat x_{0.5},\ \hat x_{0.75},\ x_{\max}$: a box from $\hat x_{0.25}$ to $\hat x_{0.75}$ with a line at the median, and whiskers out to the min and the max.

If the distribution is symmetric, the median sits in the middle of the box and the two whiskers have about the same length. Here the right whisker is much longer and the mean ($2.17$) exceeds the median ($1.6$): the sample is right-skewed, which, together with the values being nonnegative, already rules out the normal. (Some software stops the whiskers at $1.5$ box-lengths and plots the points beyond them individually, as outliers.)

5.5. Step II: Estimate the Parameters

Maximum Likelihood Estimator (MLE)

Estimator: a statistic, i.e. a function of the random observations, whose value serves to estimate some unknown parameter of interest. If the parameter is $\theta$, then the estimator is denoted by $\hat\theta$.

Ex. estimate $\mu$ (and $\sigma^2$) from $X_1, X_2, \dots, X_n$:

$$\hat\mu = \frac{X_1 + \cdots + X_n}{n} \approx \mu, \qquad S^2 = \frac{\sum_{i=1}^n (X_i - \hat\mu)^2}{n-1} \approx \sigma^2 .$$

The estimator is a function (of random variables, so it is itself a random variable). The estimate is the value of the function: once the observed values $x_1, \dots, x_n$ are plugged in, it becomes a number.

1. Construct the likelihood function $L(\theta)$, which represents the probability (likelihood) that the sequence of observations $x_1, x_2, \dots, x_n$ came from the distribution with parameter $\theta$. The data are fixed and $L$ is a function of $\theta$. Since the $X_i$’s are i.i.d. (§5.1), the joint probability is a product:
   1. Discrete case: $L(\theta) = P_\theta(X_1 = x_1, \dots, X_n = x_n) = \prod_{i=1}^n p_\theta(x_i)$
   2. Continuous case: $L(\theta) = f_\theta(x_1, \dots, x_n) = \prod_{i=1}^n f_\theta(x_i)$
2. Choose $\theta$ to maximize $L(\theta)$. The MLE is

$$\hat\theta = \arg\max_{\theta \in \Theta} L(\theta) = \arg\max_{\theta \in \Theta} \ln L(\theta), \qquad \Theta = \text{parameter set}.$$

Let $l(\theta) = \ln L(\theta)$ → the log-likelihood function. Taking $\ln$ does not move the maximum ($\ln$ is increasing), and it turns the product into a sum, which is much easier to differentiate. In smooth cases, solve $\dfrac{dl(\theta)}{d\theta} = 0$ and check $\dfrac{d^2 l(\theta)}{d\theta^2} < 0$.

Ex: $\text{Exp}(\theta)$

$$f(x) = \frac{1}{\theta}e^{-x/\theta}, \quad x \ge 0, \qquad \mathbf{E}[X] = \theta$$

(the exponential parameterized by its mean: $\theta = \beta = 1/\lambda$ of §5.2). Given $x_1, x_2, \dots, x_n$:

$$L(\theta) = \prod_{i=1}^n \frac{1}{\theta}e^{-x_i/\theta} = \frac{1}{\theta^n}\exp\Big(-\frac{1}{\theta}\sum_{i=1}^n x_i\Big),$$

$$l(\theta) = \ln L(\theta) = -n\ln\theta - \frac{1}{\theta}\sum_{i=1}^n x_i,$$

$$\frac{dl(\theta)}{d\theta} = -\frac{n}{\theta} + \frac{1}{\theta^2}\sum_{i=1}^n x_i = 0 \quad\Rightarrow\quad \hat\theta = \frac{1}{n}\sum_{i=1}^n x_i = \bar x .$$

(Check: $\dfrac{d^2 l}{d\theta^2} = \dfrac{n}{\theta^2} - \dfrac{2}{\theta^3}\sum x_i$, which at $\theta = \bar x$ equals $-\dfrac{n}{\bar x^2} < 0$: a maximum.)

The meaning: $\theta$ is the mean of the distribution, and its MLE is the sample mean, the most natural estimate.

Ex: $\text{Geo}(\theta)$ (modified)

$$p(x) = (1-\theta)^x\,\theta, \qquad x = 0, 1, 2, \dots$$

$x$ = number of trials before the first success, i.e. the number of failures (each trial is a success w.p. $\theta$). The geometric of §4.4 counts the trials up to and including the first success, $P(X = k) = (1-\theta)^{k-1}\theta$, $k = 1, 2, \dots$; here the exponent is $x$, not $x - 1$, because the modified $X$ = (usual $X$) $- 1$. $\mathbf{E}[X] = \dfrac{1-\theta}{\theta}$.

$$L(\theta) = \prod_{i=1}^n (1-\theta)^{x_i}\,\theta = \theta^n (1-\theta)^{\sum_{i=1}^n x_i},$$

$$l(\theta) = \ln L(\theta) = n\ln\theta + \Big(\sum_{i=1}^n x_i\Big)\ln(1-\theta),$$

$$\frac{dl(\theta)}{d\theta} = \frac{n}{\theta} - \frac{\sum_{i=1}^n x_i}{1-\theta} = 0 \;\Rightarrow\; n(1-\theta) = \theta\sum_{i=1}^n x_i \;\Rightarrow\; \hat\theta = \frac{n}{n + \sum_{i=1}^n x_i} = \frac{1}{1 + \bar x} .$$

($\dfrac{d^2 l}{d\theta^2} = -\dfrac{n}{\theta^2} - \dfrac{\sum x_i}{(1-\theta)^2} < 0$, so it is a maximum.)

The meaning: on average there are $\bar x$ failures before each success, i.e. one success in every $1 + \bar x$ trials, so the estimated success probability is $1/(1 + \bar x)$. Equivalently, $\mathbf{E}[X] = \frac{1-\theta}{\theta} \iff \theta = \frac{1}{1 + \mathbf{E}[X]}$, and the MLE simply replaces $\mathbf{E}[X]$ by $\bar x$.

Ex: uniform distribution $U(a, b)$

Given $x_1, x_2, \dots, x_n$; the parameter is $\theta = (a, b)$:

$$L(a, b) = \prod_{i=1}^n \frac{1}{b-a} = \frac{1}{(b-a)^n} \qquad \text{if } a \le x_i \le b \text{ for all } i,$$

and $L(a, b) = 0$ otherwise (if some $x_i$ lies outside $[a, b]$, its density is $0$). “All $x_i \in [a, b]$” means $a \le x_{(1)}$ and $b \ge x_{(n)}$.

Here setting the derivative to $0$ does not work: $1/(b-a)^n$ has no stationary point, it just keeps increasing as $b - a$ shrinks. So make $b - a$ as small as the constraint allows, i.e. push $a$ up to $x_{(1)}$ and $b$ down to $x_{(n)}$:

$$\hat a = x_{(1)} = \min_i x_i, \qquad \hat b = x_{(n)} = \max_i x_i .$$

Properties of MLEs

1. Unique for many distribution families.
2. Asymptotically unbiased: $\lim_{n \to \infty}\mathbf{E}[\hat\theta] = \theta$. For finite $n$ an MLE may be biased. E.g. for $U(a, b)$, $\hat b = X_{(n)} \le b$ always, and $\mathbf{E}[X_{(n)}] = a + (b - a)\frac{n}{n+1} < b$, but it tends to $b$ as $n \to \infty$.
3. Invariant: if instead of estimating $\theta$ we want to estimate $\phi = h(\theta)$, then the MLE is $\hat\phi = h(\hat\theta)$. E.g. the rate of the exponential is $\lambda = 1/\theta$, so $\hat\lambda = 1/\bar x$, without redoing the maximization.
4. Strongly consistent: $\lim_{n \to \infty}\hat\theta = \theta$ w.p. 1 (as in the SLLN, §2.4).
5. MLEs are asymptotically normal:

$$\sqrt n\,(\hat\theta - \theta) \xrightarrow{d} N\big(0, \delta(\theta)\big), \qquad \text{where } \delta(\theta) = \frac{-n}{\mathbf{E}\Big[\dfrac{d^2 l(\theta)}{d\theta^2}\Big]}, \quad l(\theta) = \ln L(\theta).$$

So

$$\frac{\hat\theta - \theta}{\sqrt{\delta(\theta)/n}} \xrightarrow{d} Z \sim N(0, 1), \qquad P\left(-z_{\alpha/2} \le Z \le z_{\alpha/2}\right) = 1 - \alpha .$$

A $100(1-\alpha)\%$ CI for estimating $\theta$ is

$$\left[\hat\theta - z_{\alpha/2}\sqrt{\delta(\hat\theta)/n},\;\; \hat\theta + z_{\alpha/2}\sqrt{\delta(\hat\theta)/n}\,\right].$$

How to read property 5, step by step:

• What it says. For large $n$, $\hat\theta$ is approximately normal, centred at the true $\theta$, with variance $\delta(\theta)/n$: $\hat\theta \approx N\big(\theta, \delta(\theta)/n\big)$. Standardizing (subtract the mean, divide by the standard deviation) gives the $Z$ above.
• Same steps as CI Case 1 (§2.6.1), with $\bar X(n)$ replaced by $\hat\theta$ and $\sigma^2/n$ by $\delta(\theta)/n$. Rearranging $-z_{\alpha/2} \le \frac{\hat\theta - \theta}{\sqrt{\delta(\theta)/n}} \le z_{\alpha/2}$ gives $\hat\theta - z_{\alpha/2}\sqrt{\delta(\theta)/n} \le \theta \le \hat\theta + z_{\alpha/2}\sqrt{\delta(\theta)/n}$.
• Why $\delta(\hat\theta)$ in the CI. $\delta(\theta)$ depends on the unknown $\theta$, so we plug in the estimate $\hat\theta$. For large $n$ this is fine because $\hat\theta \to \theta$ (property 4).
• Where $\delta$ comes from. $\delta(\theta)/n = -1/\mathbf{E}[l''(\theta)]$, and $l''(\theta)$ is the curvature of the log-likelihood. A sharply peaked $l$ (large $|l''|$) means the data strongly favour values near $\hat\theta$: small variance, short CI. A flat $l$ means many values of $\theta$ explain the data almost equally well: large variance, wide CI. Since $l$ is a sum of $n$ terms, $\mathbf{E}[l'']$ is proportional to $n$, so $\delta(\theta)$ does not depend on $n$, and the CI shrinks like $1/\sqrt n$.
• Check with the exponential. From the example above, $l''(\theta) = \frac{n}{\theta^2} - \frac{2}{\theta^3}\sum x_i$, and $\mathbf{E}[X_i] = \theta$, so $\mathbf{E}[l''(\theta)] = \frac{n}{\theta^2} - \frac{2n\theta}{\theta^3} = -\frac{n}{\theta^2}$ and $\delta(\theta) = \theta^2 = \operatorname{Var}(X)$. Since $\hat\theta = \bar X$, property 5 is just the CLT for $\bar X$ (§2.4), and the CI is $\bar x \pm z_{\alpha/2}\,\bar x/\sqrt n$: Case 1 of §2.6.1 with $\sigma$ replaced by its MLE $\hat\theta = \bar x$.

The log-likelihood of the exponential example with $\bar x = 2$, for $n = 10$ and $n = 50$, and the 90% CIs $\bar x \pm 1.645\,\bar x/\sqrt n$: the sharper the peak, the shorter the interval.

Ex: construct a 90% CI for $\theta$ in the modified geometric distribution $\text{Geo}(\theta)$.

$p(x) = P(X = x) = \theta(1-\theta)^x$, $x = 0, 1, 2, \dots$, $\mathbf{E}[X] = \dfrac{1-\theta}{\theta}$.

Given $x_1, \dots, x_n$ (from the example above):

$$L(\theta) = \theta^n (1-\theta)^{\sum_{i=1}^n x_i}, \qquad l(\theta) = n\ln\theta + \Big(\sum_{i=1}^n x_i\Big)\ln(1-\theta),$$

$$\frac{dl(\theta)}{d\theta} = \frac{n}{\theta} - \frac{\sum x_i}{1-\theta}, \qquad \frac{d^2 l(\theta)}{d\theta^2} = -\frac{n}{\theta^2} - \frac{\sum x_i}{(1-\theta)^2} .$$

Take expectations, using $\mathbf{E}\big[\sum X_i\big] = n\,\mathbf{E}[X] = n\frac{1-\theta}{\theta}$:

$$\mathbf{E}\Big[\frac{d^2 l(\theta)}{d\theta^2}\Big] = -\frac{n}{\theta^2} - \frac{n(1-\theta)/\theta}{(1-\theta)^2} = -\frac{n}{\theta^2} - \frac{n}{\theta(1-\theta)} = -\frac{n\big[(1-\theta) + \theta\big]}{\theta^2(1-\theta)} = -\frac{n}{\theta^2(1-\theta)} .$$

So $\delta(\theta) = \dfrac{-n}{\mathbf{E}[d^2 l/d\theta^2]} = \theta^2(1-\theta)$.

CI: with $\hat\theta = \dfrac{1}{1 + \bar x}$ and $1 - \alpha = 0.9$, $z_{\alpha/2} = z_{0.05} = 1.645$:

$$\left[\hat\theta - 1.645\sqrt{\frac{\hat\theta^2(1-\hat\theta)}{n}},\;\; \hat\theta + 1.645\sqrt{\frac{\hat\theta^2(1-\hat\theta)}{n}}\,\right].$$

E.g. $n = 50$ observations with $\bar x = 2$: $\hat\theta = 1/3$, $\delta(\hat\theta) = \frac19 \cdot \frac23 = 0.0741$, $\sqrt{\delta(\hat\theta)/n} = 0.0385$, so the 90% CI is $0.333 \pm 1.645 \times 0.0385 = [0.270,\ 0.397]$.

5.6. Step III: Determine the Quality of Fit

1. Graphical methods
   1. quantile–quantile plot (Q–Q plot)
   2. probability–probability plot (P–P plot)
2. Statistical methods
   1. chi-square test
   2. Kolmogorov–Smirnov (K–S) test

Running example (used for the plots and the test below). $n = 20$ interarrival times, sorted:

$$0.16,\ 0.18,\ 0.19,\ 0.38,\ 0.41,\ 0.90,\ 1.05,\ 1.08,\ 1.19,\ 1.22,\ 1.32,\ 1.38,\ 1.43,\ 1.97,\ 2.09,\ 2.19,\ 2.29,\ 2.53,\ 2.55,\ 2.95 .$$

Step I: hypothesize an exponential distribution. Step II: the MLE of its mean is $\hat\theta = \bar x = 27.46/20 = 1.373$, so $\hat F(x) = 1 - e^{-x/1.373}$. Step III asks how well this $\hat F$ fits. (The numbers were simulated from an exponential distribution, so the exponential family is in fact correct here.)

Graphical methods

Q–Q plot

Let $\hat F$ be the fitted theoretical distribution and $\tilde F$ the empirical distribution based on $x_1, \dots, x_n$. A Q–Q plot is a graph of the $q_i$-quantile of $\hat F$ vs the $q_i$-quantile of $\tilde F$.

$$F(x) = \frac{\text{number of } x_i\text{'s} \le x}{n}, \qquad \tilde F(x) = F(x) - \frac{0.5}{n}, \qquad q_i = \frac{i - 0.5}{n}, \quad i = 1, \dots, n .$$

(This $\tilde F$, the empirical c.d.f. shifted down by $0.5/n$, is not the interpolated $\tilde F$ of §5.3.)

$\tilde F(x_{(1)}) = \frac{1}{n} - \frac{0.5}{n} = \frac{1 - 0.5}{n} = q_1$, and in general $\tilde F(x_{(i)}) = \frac{i - 0.5}{n} = q_i$ ⇒ the $q_i$-quantile of $\tilde F$ is just $x_{(i)}$.

Why $0.5/n$? Without this term, $q_n = n/n = 1$, and we would need the 1-quantile of $\hat F$. For a distribution with infinite range, like the normal or the exponential, $\hat F^{-1}(1) = +\infty$: not possible. With the $0.5/n$, every $q_i$ lies strictly inside $(0, 1)$, and symmetrically: $q_1 = \frac{0.5}{n}$, $q_n = 1 - \frac{0.5}{n}$.

Plot $x^M_{q_i} = \hat F^{-1}(q_i)$ (model quantile) as a function of $x^S_{q_i} = x_{(i)}$ (sample quantile), $i = 1, \dots, n$. If $\hat F$ fits well, the points lie close to the $45°$ line $y = x$.

Left: at height $q_i$, read the sample quantile $x_{(i)}$ off the staircase $\tilde F$ and the model quantile $x^M = \hat F^{-1}(q_i)$ off the curve $\hat F$. The Q–Q plot compares these two x-values (a horizontal gap). Right: the Q–Q plot of the running example, with $\hat F^{-1}(q) = -1.373\ln(1-q)$. The points follow the line except at the upper end: the two largest observations, $2.55$ and $2.95$, are much smaller than the model quantiles $3.56$ and $5.06$, i.e. the sample has a shorter right tail than the fitted exponential.

Used for continuous distributions (it needs the quantile function $\hat F^{-1}$; for a discrete $\hat F$ many $q_i$ share the same quantile, and the plot shows little).

P–P plot

A P–P plot is a graph of the model probability $\hat F(x_{(i)})$ vs the sample probability $\tilde F(x_{(i)}) = q_i = \frac{i - 0.5}{n}$, for $i = 1, \dots, n$. Used for both continuous and discrete data sets. If $\hat F$ fits well, the points lie close to the line $y = x$ in $[0, 1]^2$.

Left: at $x_{(i)}$, read the probability on the curve $\hat F$ and on the staircase $\tilde F$. The P–P plot compares these two probabilities (a vertical gap). Right: the P–P plot of the running example. The largest deviations are in the middle, e.g. $\hat F(x_{(6)}) = \hat F(0.90) = 0.48$ vs $q_6 = 0.275$ (the sample has no observations between $0.41$ and $0.90$). At the upper end the points are much closer to the line: the short right tail that stands out in the Q–Q plot is barely visible here.

Q–Q or P–P?

A small change in the tail causes a large gap in the Q–Q plot. In the tails $\hat F$ is flat (the density is small), so a tiny difference in probability corresponds to a large difference in $x$; the Q–Q plot compares x-values, so differences in the tails are magnified. In the middle $\hat F$ is steep, so a small difference in $x$ corresponds to a large difference in probability; the P–P plot compares probabilities, so differences in the middle are magnified. (In the tails both $\hat F$ and $\tilde F$ are close to $0$ or $1$, so the P–P points crowd into the corners and tail differences are hidden.)

• If you want to know whether $\hat F$ provides a good approximation of $\tilde F$ in the tails → Q–Q plot.
• If you want to check the middle (the bulk) of the distribution → P–P plot. It is also the one of the two that works for discrete data.

The running example shows both effects on the same data: the Q–Q plot’s biggest gap is in the right tail, the P–P plot’s in the middle.

Goodness-of-fit tests (hypothesis tests)

They test the null hypothesis

$$H_0: \text{the } X_i\text{'s are i.i.d. with distribution } \hat F .$$

As with the RNG tests (§3.4), “do not reject $H_0$” does not prove that $\hat F$ is correct, only that the data show no evidence against it. The sample size matters: with little data almost nothing is rejected, and with a lot of data almost every $\hat F$ is rejected, since no theoretical distribution fits real data exactly.

Chi-square test

1. Split the range of $\hat F$ (not the range of $x_1, \dots, x_n$) into $k$ adjacent intervals $[a_0, a_1), [a_1, a_2), \dots, [a_{k-1}, a_k)$. It is possible that $a_0 = -\infty$ and/or $a_k = +\infty$. (Every value that $\hat F$ can produce must fall in some interval, so that the $p_j$ in step 3 add up to 1.)
2. Count the number of $x_i$’s in each interval: $f_j$ = number of $x_i$’s in the $j$-th interval, $\sum_{j=1}^k f_j = n$.
3. Assume $H_0$ is true; compute the expected proportion $p_j$ of $x_i$’s that would fall in the $j$-th interval:

$$p_j = \int_{a_{j-1}}^{a_j} \hat f(x)\,dx = \hat F(a_j) - \hat F(a_{j-1}) \quad \text{(continuous)}, \qquad p_j = \sum_{a_{j-1} \le x < a_j} \hat p(x) \quad \text{(discrete)},$$

where in the discrete case the sum runs over the possible values $x$ of the fitted distribution (not over the data). The expected number in interval $j$ is $np_j$.
4. Test statistic

$$\chi^2 = \sum_{j=1}^k \frac{(f_j - np_j)^2}{np_j} .$$

Reject $H_0$ at significance level $\alpha$ if $\chi^2 > \chi^2_{k-1,\,\alpha}$ (see the $\chi^2$ table).

Remarks.

• It is the same statistic as the RNG uniformity test of §3.4.1, which is the special case $\hat F = U(0,1)$ with $k$ equal intervals, $p_j = 1/k$.
• Choosing the intervals: a common choice is equiprobable intervals, $p_j = 1/k$, i.e. $a_j = \hat F^{-1}(j/k)$, with $np_j \ge 5$ (as in §3.4.1).
• Degrees of freedom. If $m$ parameters of $\hat F$ were estimated from the same data (e.g. by MLE), the distribution of $\chi^2$ under $H_0$ lies between $\chi^2_{k-m-1}$ and $\chi^2_{k-1}$. Using $\chi^2_{k-1,\alpha}$, as in step 4, is conservative: the actual probability of a type I error is at most $\alpha$. So if $\chi^2 \le \chi^2_{k-m-1,\alpha}$, do not reject; if $\chi^2 > \chi^2_{k-1,\alpha}$, reject; in between, the decision depends on which critical value is used.

Example (running example, $\hat F$ = exponential with $\hat\theta = 1.373$). Take $k = 4$ equiprobable intervals, so $np_j = 20/4 = 5$, with boundaries $a_j = \hat F^{-1}(j/4) = -1.373\ln(1 - j/4)$:

$j$	interval	$f_j$	$np_j$	$(f_j - np_j)^2/np_j$
1	$[0,\ 0.395)$	4	5	0.2
2	$[0.395,\ 0.952)$	2	5	1.8
3	$[0.952,\ 1.903)$	7	5	0.8
4	$[1.903,\ \infty)$	7	5	0.8
	total	20	20	$\chi^2 = 3.6$

$m = 1$ parameter was estimated, so the critical value at $\alpha = 0.05$ is between $\chi^2_{2,\,0.05} = 5.991$ and $\chi^2_{3,\,0.05} = 7.815$. Since $\chi^2 = 3.6 < 5.991$, do not reject $H_0$ either way: no evidence against the exponential fit.

Note what the test cannot see: interval 4 only counts how many observations exceed $1.903$ (7, close to the expected 5), not how far out they go, so the short right tail visible in the Q–Q plot does not affect $\chi^2$. With $n = 20$ that short tail is a chance feature of the sample anyway, since the data do come from an exponential.