2. Review of Probability and Statistics

A recap of the probability and statistics needed for simulation: probability axioms, random variables, convergence modes, stochastic processes, and estimation of means and variances.

2.1. Axioms of Probability

• Sample space $\Omega$: the collection of all outcomes of a random experiment. Each element $\omega \in \Omega$ is called a sample point.
• Event: any subset of the sample space, $A \subseteq \Omega$.
• Probability: a function that assigns to each event $A$ a number $P(A)$. It is not an arbitrary function; it must satisfy the following axioms:
  1. $0 \le P(A) \le 1$
  2. $P(\Omega) = 1$
  3. If $A_1, A_2, \dots$ are pairwise disjoint ($A_i \cap A_j = \emptyset$ for $i \ne j$), then $P\left(\bigcup_{i=1}^\infty A_i\right) = \sum_{i=1}^\infty P(A_i)$.

2.2. Random Variables and Their Properties

A random variable is a function that assigns to each $\omega \in \Omega$ a numerical value: $X : \Omega \to (-\infty, \infty)$.

The cumulative distribution function (c.d.f.) is

$$F(x) = P(X \le x) = P(\{\omega : X(\omega) \le x\}).$$

Properties:

1. $F$ is non-decreasing in $x$, and right-continuous with left limits. At a jump, the value at the point equals the value after the jump: $F(x_0) = F(x_0^+)$.

2. $\lim_{x \to \infty} F(x) = 1$, $\lim_{x \to -\infty} F(x) = 0$

3. $P(a < X \le b) = F(b) - F(a)$

2.2.1. Two types of random variables

Discrete: $X$ takes a countable number of different values. A discrete r.v. is described by its probability mass function (p.m.f.) $p_X(x) = P(X = x)$:

1. $0 \le p_X(x) \le 1$
2. $\sum_{\text{all } x} p_X(x) = 1$
3. $F(x) = P(X \le x) = \sum_{y \le x} p_X(y)$
4. $p_X(x) = F(x) - F(x^-)$ (value right after the jump minus value right before the jump)

Continuous: there exists a function $f$ (the probability density function, p.d.f.) such that $F(x) = \int_{-\infty}^{x} f(y)\,dy$:

1. $f(x) \ge 0$
2. $\int_{-\infty}^{\infty} f(x)\,dx = 1$
3. $P(a \le X \le b) = \int_a^b f(x)\,dx$
4. $\dfrac{dF(x)}{dx} = f(x)$

2.2.2. Expectation and variance

Expectation of a r.v.:

$$\mathbf{E}[X] = \begin{cases} \displaystyle\sum_{\text{all } x} x\, p_X(x) & X \text{ discrete} \\[6pt] \displaystyle\int_{-\infty}^{\infty} x\, f(x)\,dx & X \text{ continuous} \end{cases}$$

Expectation of a function of a r.v., $Y = g(X)$:

$$\mathbf{E}[Y] = \begin{cases} \displaystyle\sum_{\text{all } y} y\, p_Y(y) = \sum_{\text{all } x} g(x)\, p_X(x) & X \text{ discrete} \\[6pt] \displaystyle\int_{-\infty}^{\infty} y\, f_Y(y)\,dy = \int_{-\infty}^{\infty} g(x)\, f_X(x)\,dx & X \text{ continuous} \end{cases}$$

Variance:

$$\begin{aligned} \operatorname{Var}(X) &= \mathbf{E}\left[(X - \mathbf{E}[X])^2\right] \\ &= \begin{cases} \displaystyle\sum_{\text{all } x} (x - \mathbf{E}[X])^2\, p_X(x) & X \text{ discrete} \\[6pt] \displaystyle\int_{-\infty}^{\infty} (x - \mathbf{E}[X])^2\, f(x)\,dx & X \text{ continuous} \end{cases} \\ &= \mathbf{E}[X^2] - \left(\mathbf{E}[X]\right)^2 \end{aligned}$$

Standard deviation: $\sqrt{\operatorname{Var}(X)}$.

2.3. Bivariate Random Variables

Let $X : \Omega \to \mathbb{R}$ and $Y : \Omega \to \mathbb{R}$ be two r.v.s. The pair $(X, Y)$ is called a bivariate r.v., described by the joint c.d.f.

$$F(x, y) = P(X \le x,\; Y \le y).$$

Marginal c.d.f.s

$$\begin{aligned} F_X(x) &= P(X \le x) = \lim_{y \to \infty} F(x, y), \\ F_Y(y) &= P(Y \le y) = \lim_{x \to \infty} F(x, y). \end{aligned}$$

Discrete case. Joint p.m.f. and marginals:

$$\begin{aligned} p(x, y) &= P(X = x,\; Y = y), \\ p_X(x) &= \sum_{\text{all } y} p(x, y), \qquad p_Y(y) = \sum_{\text{all } x} p(x, y). \end{aligned}$$

Conditional p.m.f. of $X$ given $Y = y$:

$$p_{X \mid Y}(x \mid y) = \frac{p(x, y)}{p_Y(y)}, \qquad p_Y(y) > 0.$$

Continuous case. Joint p.d.f. $f$ and marginals:

$$\begin{aligned} F(x, y) &= \int_{-\infty}^{x} \int_{-\infty}^{y} f(u, v)\,dv\,du, \\ f_X(x) &= \int_{-\infty}^{\infty} f(x, y)\,dy, \qquad f_Y(y) = \int_{-\infty}^{\infty} f(x, y)\,dx. \end{aligned}$$

Conditional p.d.f. of $X$ given $Y = y$:

$$f_{X \mid Y}(x \mid y) = \frac{f(x, y)}{f_Y(y)}, \qquad f_Y(y) > 0.$$

Independence. $X$ and $Y$ are independent iff $F(x, y) = F_X(x)\, F_Y(y)$ for all $x, y$; equivalently

• discrete: $p(x, y) = p_X(x)\, p_Y(y)$
• continuous: $f(x, y) = f_X(x)\, f_Y(y)$

Expectation of $Z = g(X, Y)$:

$$\mathbf{E}[Z] = \mathbf{E}[g(X, Y)] = \begin{cases} \displaystyle\sum_{\text{all } x}\sum_{\text{all } y} g(x, y)\, p(x, y) & \text{discrete} \\[6pt] \displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} g(x, y)\, f(x, y)\,dx\,dy & \text{continuous} \end{cases}$$

Covariance

$$\begin{aligned} \operatorname{cov}(X, Y) &= \mathbf{E}\left[(X - \mathbf{E}[X])(Y - \mathbf{E}[Y])\right] \\ &= \mathbf{E}[XY] - \mathbf{E}[X]\,\mathbf{E}[Y]. \end{aligned}$$

• If $X = Y$, then $\operatorname{cov}(X, Y) = \operatorname{Var}(X)$.
• If $X$ and $Y$ are independent, then $\operatorname{cov}(X, Y) = 0$, i.e. $\mathbf{E}[XY] = \mathbf{E}[X]\,\mathbf{E}[Y]$.
• The converse fails in general: $\operatorname{cov}(X, Y) = 0 \not\Rightarrow$ independence. But if $(X, Y)$ is jointly normal, then $\operatorname{cov}(X, Y) = 0 \iff X$ and $Y$ are independent.
• $\operatorname{cov} > 0$: positively correlated — when $X$ is above its mean, $Y$ tends to be above its mean (and below with below).
• $\operatorname{cov} < 0$: negatively correlated — when $X$ is above its mean, $Y$ tends to be below its mean.

Correlation coefficient

$$\operatorname{corr}(X, Y) = \frac{\operatorname{cov}(X, Y)}{\sqrt{\operatorname{Var}(X)}\,\sqrt{\operatorname{Var}(Y)}} \in [-1, 1].$$

2.4. Convergence Modes

Let $X_1, X_2, \dots$ be a sequence of r.v.s. Each $X_n$ is a function $X_n(\omega)$, so “converge to $X^*(\omega)$” is not the ordinary limit of a sequence of numbers like $\lim_{n \to \infty} 1/n = 0$; we need to say in what sense the functions converge.

Convergence with probability 1 (a.s. convergence, strong convergence).
Here $X_n$ is still a random variable, so the statement is a probability:

$$P\left(\left\{\omega : \lim_{n \to \infty} X_n(\omega) = X^*(\omega)\right\}\right) = 1.$$

For almost every sampled $\omega$, plugging it into $X_n$ gives a sequence of numbers that converges to $X^*(\omega)$. The equality may fail on a set of $\omega$ of probability zero (e.g. on $\Omega = [0,1]$ with the uniform measure, failing only at $\omega = 0.5$ still gives probability 1).

Convergence in distribution (weak convergence).
Let $F_n$ be the c.d.f. of $X_n$ and $F^*$ the c.d.f. of $X^*$. These are ordinary functions, not random, so we take the limit directly:

$$\lim_{n \to \infty} F_n(x) = F^*(x) \quad \text{at all continuity points } x \text{ of } F^*.$$

Theorem (Strong Law of Large Numbers, SLLN).
Let $X_1, X_2, \dots$ be i.i.d. with $\mathbf{E}|X_1| < \infty$ and $\mathbf{E}[X_1] = \mu$. Let $\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i$. Then $\bar{X}_n \to \mu$ with probability 1:

$$P\left(\left\{\omega : \lim_{n \to \infty} \frac{1}{n}\sum_{i=1}^{n} X_i(\omega) = \mu\right\}\right) = 1.$$

This is exactly the first mode, with limit $X^* \equiv \mu$ a constant random variable. The only assumption is that the mean exists; the variance may be infinite.

(Contrast the Weak LLN, which only gives convergence in probability: $P(|\bar{X}_n - \mu| > \varepsilon) \to 0$ for each fixed $\varepsilon > 0$. This is weaker — it does not say that any individual sample path settles down.)

Theorem (Central Limit Theorem, CLT).
Let $X_1, X_2, \dots$ be i.i.d. with $\mathbf{E}[X_1] = \mu$ and $0 < \operatorname{Var}(X_1) = \sigma^2 < \infty$. Define

$$Z_n = \frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}}, \qquad F_n(x) = P(Z_n \le x).$$

Then $Z_n \to N(0, 1)$ in distribution:

$$\lim_{n \to \infty} F_n(x) = \Phi(x) \quad \text{for all } x,$$

where $\Phi$ is the $N(0,1)$ c.d.f. (Since $\Phi$ is continuous everywhere, every $x$ is a continuity point.)

Why these two are the canonical examples. Same setup, two different questions, two different modes:

	SLLN	CLT
object	$\bar{X}_n$ itself	rescaled fluctuation $\sqrt{n}(\bar{X}_n - \mu)/\sigma$
mode	with probability 1 (strong)	in distribution (weak)
limit	constant $\mu$	$N(0, 1)$
needs	$\mathbf{E}\lvert X_1\rvert < \infty$	$\operatorname{Var}(X_1) < \infty$

SLLN says the average lands on $\mu$; CLT says how it wobbles on the way in — at scale $1/\sqrt{n}$. CLT cannot be strengthened to a.s. convergence: for any fixed $\omega$, $Z_n(\omega)$ does not settle down but keeps oscillating forever. Only its distribution stabilizes. This illustrates that weak convergence is genuinely weaker, not just a technical relaxation.

2.5. Stochastic Processes

A stochastic process is a collection of r.v.s $\{X(t) : t \in T\}$ indexed by $t$, which usually represents time.

• If $T$ is discrete / continuous, we have a discrete-time / continuous-time process.
• If the values $X(t)$ take are discrete / continuous, we have a discrete-state / continuous-state process.

Example: single-server FCFS queue. Interarrival times $X_1, X_2, \dots$ are i.i.d. and service times $Y_1, Y_2, \dots$ are i.i.d. Let $D_i$ be the delay of the $i$-th customer. Then $D_1 = 0$ and

$$D_{n+1} = \max\{D_n + Y_n - X_{n+1},\; 0\} \qquad \text{(Lindley's equation)}.$$

For instance $D_2 = \max\{Y_1 - X_2, 0\}$: customer 2 waits only if customer 1’s service outlasts the gap between the two arrivals. $\{D_n\}$ is a discrete-time, continuous-state process.

Example: Poisson process. Let $S_n$ be the arrival time of the $n$-th customer, $S_0 = 0$, and let the interarrival times $X_n = S_n - S_{n-1}$ be i.i.d. $\text{Exp}(\lambda)$, so $S_n = \sum_{i=1}^n X_i$. The number of arrivals up to time $t$ is

$$N(t) = \max\{n : S_n \le t\}.$$

$\{N(t) : t \ge 0\}$ is a continuous-time, discrete-state process, called a Poisson process with rate $\lambda$.

2.6. Estimation of Means and Variances

Let $X_1, \dots, X_n$ be i.i.d. with mean $\mu$ and variance $\sigma^2$.

Sample mean

$$\bar{X}(n) = \frac{X_1 + X_2 + \cdots + X_n}{n}, \qquad \mathbf{E}[\bar{X}(n)] = \mu,$$

so $\bar{X}(n)$ is an unbiased estimator of $\mu$.

Sample variance

$$S^2(n) = \frac{1}{n-1}\sum_{i=1}^n \left(X_i - \bar{X}(n)\right)^2, \qquad \mathbf{E}[S^2(n)] = \sigma^2,$$

so $S^2(n)$ is an unbiased estimator of $\sigma^2$.

Both $\bar{X}(n)$ and $S^2(n)$ are point estimators, and both are themselves random variables.

Idea. Replace the point estimator by an interval estimator (confidence interval), so that one can assert with a certain degree of confidence that the interval contains the true mean $\mu$.

2.6.1. Confidence Intervals

Case 1. Population distribution arbitrary, $\sigma^2$ known.

Notation. Let $Z \sim N(0,1)$ and define the upper $\alpha$-quantile $z_\alpha$ by

$$P(Z > z_\alpha) = \alpha.$$

By symmetry of the standard normal, $P(Z < -z_\alpha) = \alpha$, and hence

$$P\left(-z_{\alpha/2} \le Z \le z_{\alpha/2}\right) = 1 - \alpha.$$

Derivation. By the Central Limit Theorem, for large $n$,

$$\frac{\bar{X}(n) - \mu}{\sigma / \sqrt{n}} \;\approx\; Z \sim N(0,1).$$

Therefore

$$P\left(-z_{\alpha/2} \le \frac{\bar{X}(n) - \mu}{\sigma/\sqrt{n}} \le z_{\alpha/2}\right) \;\approx\; P\left(-z_{\alpha/2} \le Z \le z_{\alpha/2}\right) = 1 - \alpha.$$

Rearranging the event inside the probability,

$$-z_{\alpha/2} \le \frac{\bar{X}(n) - \mu}{\sigma/\sqrt{n}} \le z_{\alpha/2} \iff \bar{X}(n) - z_{\alpha/2}\frac{\sigma}{\sqrt{n}} \le \mu \le \bar{X}(n) + z_{\alpha/2}\frac{\sigma}{\sqrt{n}}.$$

Conclusion. An approximate $100(1-\alpha)\%$ confidence interval for $\mu$ is

$$\left[\, \bar{X}(n) - z_{\alpha/2}\frac{\sigma}{\sqrt{n}},\; \bar{X}(n) + z_{\alpha/2}\frac{\sigma}{\sqrt{n}} \,\right].$$

Case 2. Population normal, $\sigma^2$ unknown.

Notation. Let $T \sim t_{n-1}$ (Student’s $t$ with $n-1$ degrees of freedom) and define $t_{\alpha,\,n-1}$ by

$$P(T > t_{\alpha,\,n-1}) = \alpha.$$

Since the $t$ distribution is symmetric about $0$,

$$P\left(-t_{\alpha/2,\,n-1} \le T \le t_{\alpha/2,\,n-1}\right) = 1 - \alpha.$$

Derivation. Let $X_1, \dots, X_n \overset{\text{i.i.d.}}{\sim} N(\mu, \sigma^2)$ with $\sigma^2$ unknown, and estimate $\sigma^2$ by $S^2(n)$. Replacing $\sigma$ by $S(n)$ in the standardized mean gives

$$\frac{\bar{X}(n) - \mu}{S(n)/\sqrt{n}} \;\sim\; t_{n-1}$$

exactly (not just approximately), because $\bar{X}(n)$ and $S^2(n)$ are independent under normality, with $\frac{(n-1)S^2(n)}{\sigma^2} \sim \chi^2_{n-1}$. Hence

$$P\left(-t_{\alpha/2,\,n-1} \le \frac{\bar{X}(n) - \mu}{S(n)/\sqrt{n}} \le t_{\alpha/2,\,n-1}\right) = 1 - \alpha,$$

and rearranging as in Case 1 gives:

Conclusion. An exact $100(1-\alpha)\%$ confidence interval for $\mu$ is

$$\left[\, \bar{X}(n) - t_{\alpha/2,\,n-1}\frac{S(n)}{\sqrt{n}},\; \bar{X}(n) + t_{\alpha/2,\,n-1}\frac{S(n)}{\sqrt{n}} \,\right].$$

Comparison of the two intervals.

For every $n$ and every $\alpha \in (0,1)$,

$$t_{\alpha/2,\,n-1} \;>\; z_{\alpha/2},$$

because the $t_{n-1}$ distribution has heavier tails than $N(0,1)$. Consequently the $t$-interval is wider than the $z$-interval with the same $S(n)$: the extra width is the price of estimating $\sigma$ from the data rather than knowing it.

As $n \to \infty$, $t_{n-1} \to N(0,1)$ and so $t_{\alpha/2,\,n-1} \downarrow z_{\alpha/2}$; the two intervals coincide in the limit. In practice the difference is negligible for $n \gtrsim 30$.

Which is “better” depends on what is known:

• $\sigma^2$ known (Case 1): use $z_{\alpha/2}$. The interval is narrower, but only approximately at level $1-\alpha$ unless the population is normal.
• $\sigma^2$ unknown and population normal (Case 2): use $t_{\alpha/2,\,n-1}$. The interval is exactly at level $1-\alpha$; substituting $z_{\alpha/2}$ here would give an interval that is too narrow and under-covers $\mu$, especially for small $n$.

2.6.2. Hypothesis Testing Concerning the Mean

Goal: determine whether a claim about $\mu$ is true or not.

$$H_0: \mu = \mu_0 \quad \text{(null hypothesis)} \qquad \text{vs.} \qquad H_1: \mu \ne \mu_0 \quad \text{(alternative hypothesis)}$$

Steps in hypothesis testing

1. Set up $H_0$ and $H_1$.
2. Specify the significance level $\alpha$: the maximum allowed type I error probability, $P(\text{reject } H_0 \mid H_0 \text{ true}) \le \alpha$.

	Accept $H_0$	Reject $H_0$
$H_0$ is true	correct	type I error
$H_1$ is true	type II error	correct

3. Choose the test statistic and the rejection region.
   For a normal population with $\sigma^2$ unknown, use

   $$t_n = \frac{\bar{X}(n) - \mu_0}{S(n)/\sqrt{n}},$$

   which has a $t_{n-1}$ distribution under $H_0$. The rejection region is $\{\,|t_n| > t_{\alpha/2,\,n-1}\,\}$, so that $P(\text{reject } H_0 \mid H_0) = \alpha$.
4. Substitute the observed values into the test statistic and make the decision:
   • if $|t_n| > t_{\alpha/2,\,n-1}$, reject $H_0$;
   • if $|t_n| \le t_{\alpha/2,\,n-1}$, do not reject $H_0$.

Remark. Rejecting $H_0$ is equivalent to $\mu_0 \notin \left[\bar{X}(n) \pm t_{\alpha/2,\,n-1}\, S(n)/\sqrt{n}\right]$, the CI from Case 2. If instead $\sigma^2$ is known (Case 1), replace $S(n)$ by $\sigma$ and $t_{\alpha/2,\,n-1}$ by $z_{\alpha/2}$.

2.6.3. Statistical tables ($z$, $t$, $\chi^2$)

Every test so far ends by comparing a statistic against a threshold read from a table. All three tables answer the same question in the same way.

What a table contains. For a distribution with c.d.f. $F$, the upper $\alpha$-quantile $x_\alpha$ is defined by

$$P(X > x_\alpha) = \alpha, \qquad\text{equivalently}\qquad F(x_\alpha) = 1 - \alpha .$$

So $x_\alpha = F^{-1}(1-\alpha)$: the table is just the inverse c.d.f. evaluated at a few convenient values of $\alpha$. This is why $z_\alpha$, $t_{\alpha, n-1}$, and $\chi^2_{d,\alpha}$ all have the same shape of definition — only the distribution changes.

How the numbers are computed. There is no closed form for $F^{-1}$ in any of these three cases, so the entries are computed numerically: evaluate $F(x)$ (itself an integral, computed by series expansion or numerical quadrature) and solve $F(x) = 1-\alpha$ for $x$ by root-finding. For the chi-square distribution with $d$ degrees of freedom the density is

$$f_d(x) = \frac{1}{2^{d/2}\,\Gamma(d/2)}\, x^{d/2 - 1} e^{-x/2}, \qquad x > 0,$$

and $\chi^2_{d,\alpha}$ is the number that leaves area $\alpha$ to its right under this curve. Before computers these values were tabulated once and printed in the back of textbooks; today any software returns them directly (in R, qchisq(1-alpha, df); in Python, scipy.stats.chi2.ppf(1-alpha, df)).

Reading the $\chi^2$ table. Rows are degrees of freedom $d$, columns are $\alpha$. The entry is $\chi^2_{d,\alpha}$.

$d \backslash \alpha$	0.10	0.05	0.025	0.01
1	2.706	3.841	5.024	6.635
2	4.605	5.991	7.378	9.210
3	6.251	7.815	9.348	11.345
4	7.779	9.488	11.143	13.277
5	9.236	11.070	12.833	15.086
6	10.645	12.592	14.449	16.812
8	13.362	15.507	17.535	20.090
10	15.987	18.307	20.483	23.209
15	22.307	24.996	27.488	30.578
20	28.412	31.410	34.170	37.566
24	33.196	36.415	39.364	42.980

E.g. $\chi^2_{3,\,0.05} = 7.815$: a $\chi^2_3$ variable exceeds $7.815$ with probability $0.05$.

Two things to notice.

• The entries grow with $d$, because $\mathbf{E}[\chi^2_d] = d$ and $\operatorname{Var}(\chi^2_d) = 2d$: more terms in the sum of squares means a larger typical value. A raw $\chi^2$ value therefore means nothing without its degrees of freedom.
• The $\chi^2$ table is one-sided (upper tail only), unlike the $z$ and $t$ tables which are used two-sided for confidence intervals. That is because $\chi^2$ statistics measure squared deviation, so only large values are evidence against $H_0$.

For large $d$ a table is unnecessary: $\chi^2_d \approx N(d, 2d)$, so $\chi^2_{d,\alpha} \approx d + z_\alpha\sqrt{2d}$. Check at $d = 24$, $\alpha = 0.05$: $24 + 1.645\sqrt{48} = 35.4$ against the exact $36.415$.