1 1. Aim and teaching plan

This lecture assumes that the probability and inference theory is already known. The goal is to translate that theory into a coherent R workflow.

The session has four connected themes:

  1. probability distributions and the d/p/q/r grammar in R;
  2. one reusable classroom dataset for z, t, chi-square, and F procedures;
  3. a visual demonstration of Type I error, Type II error, and statistical power;
  4. simulation of Brownian motion and Brownian bridge, ending with the Kolmogorov–Smirnov goodness-of-fit test.

A useful 60-minute allocation is:

Time Topic
0–5 min R distribution grammar: d, p, q, r
5–12 min Discrete: uniform, binomial, Poisson
12–20 min Continuous: uniform, normal, gamma, beta
20–25 min Create one reusable Section A/B score dataset
25–31 min z-test + confidence interval
31–36 min t distribution + t-tests + confidence intervals
36–41 min Chi-square distribution + test
41–45 min F distribution + variance-ratio test + CI
45–53 min Type I/II errors and dynamic power demonstration
53–57 min Brownian motion + Brownian bridge
57–60 min Kolmogorov distribution idea + KS goodness-of-fit

2 2. Packages

Base R already provides nearly everything needed. For a convenient z-test function we use the CRAN package BSDA, whose z.test() function handles one- and two-sample z procedures when the population standard deviation(s) are treated as known.

Run the installation command only once on a computer.

install.packages("BSDA")

Then load it in each new R session:

if (requireNamespace("BSDA", quietly = TRUE)) {
  library(BSDA)
} else {
  message("Package 'BSDA' is not installed. Run install.packages('BSDA') before the z-test section.")
}

3 3. The R distribution grammar: learn it once

For most standard probability distributions, R uses four prefixes:

Prefix Meaning Mathematical role
d density / probability mass \(f(x)\) or \(P(X=x)\)
p cumulative probability \(F(x)=P(X\le x)\)
q quantile \(F^{-1}(p)\)
r random generation simulate observations

For the normal distribution:

dnorm(0)                 # density at x = 0
#> [1] 0.3989423
pnorm(1.96)              # P(Z <= 1.96)
#> [1] 0.9750021
qnorm(0.975)             # 97.5th percentile
#> [1] 1.959964
set.seed(101)
rnorm(5)                 # five N(0,1) observations
#> [1] -0.3260365  0.5524619 -0.6749438  0.2143595  0.3107692

The same grammar appears repeatedly:

Normal:     dnorm   pnorm   qnorm   rnorm
Binomial:   dbinom  pbinom  qbinom  rbinom
Poisson:    dpois   ppois   qpois   rpois
Uniform:    dunif   punif   qunif   runif
Gamma:      dgamma  pgamma  qgamma  rgamma
Beta:       dbeta   pbeta   qbeta   rbeta
t:          dt      pt      qt      rt
Chi-square: dchisq  pchisq  qchisq  rchisq
F:          df      pf      qf      rf

4 4. Discrete distributions

4.1 4.1 Discrete uniform distribution

A fair die is a simple discrete-uniform model:

\[ P(X=x)=\frac16,\qquad x=1,\ldots,6. \]

x_die <- 1:6
pmf_die <- rep(1/6, 6)

barplot(
  pmf_die,
  names.arg = x_die,
  ylim = c(0, 0.22),
  col = "#0072B2",
  border = NA,
  xlab = "Outcome",
  ylab = "Probability",
  main = "Discrete Uniform Distribution: Fair Die"
)

set.seed(110)
die_sim <- sample(1:6, size = 1000, replace = TRUE)

barplot(
  prop.table(table(die_sim)),
  ylim = c(0, 0.22),
  col = "#56B4E9",
  border = NA,
  xlab = "Outcome",
  ylab = "Relative frequency",
  main = "1000 Simulated Die Rolls"
)
abline(h = 1/6, lty = 2, lwd = 2)

The second plot demonstrates the empirical frequencies approaching the theoretical probabilities.

4.2 4.2 Binomial distribution

If \(X\sim\mathrm{Binomial}(n,p)\),

\[ P(X=x)=\binom{n}{x}p^x(1-p)^{n-x}. \]

x_bin <- 0:20
bin_pmf <- dbinom(x_bin, size = 20, prob = 0.40)

barplot(
  bin_pmf,
  names.arg = x_bin,
  col = "#009E73",
  border = NA,
  xlab = "Number of successes",
  ylab = "Probability",
  main = "Binomial(20, 0.40)"
)

dbinom(8, size = 20, prob = 0.40)  # P(X = 8)
#> [1] 0.1797058
pbinom(8, size = 20, prob = 0.40)  # P(X <= 8)
#> [1] 0.5955987
qbinom(0.95, size = 20, prob = 0.40)
#> [1] 12
set.seed(120)
rbinom(10, size = 20, prob = 0.40)
#>  [1]  7  6  7 10  6  8 13  5 10 11

4.3 4.3 Poisson distribution

For \(X\sim\mathrm{Poisson}(\lambda)\),

\[ P(X=x)=e^{-\lambda}\frac{\lambda^x}{x!}. \]

x_pois <- 0:15
pois_pmf <- dpois(x_pois, lambda = 4)

barplot(
  pois_pmf,
  names.arg = x_pois,
  col = "#D55E00",
  border = NA,
  xlab = "Count",
  ylab = "Probability",
  main = "Poisson Distribution: lambda = 4"
)

dpois(5, lambda = 4)
#> [1] 0.1562935
ppois(5, lambda = 4)
#> [1] 0.7851304
qpois(0.95, lambda = 4)
#> [1] 8
set.seed(130)
rpois(10, lambda = 4)
#>  [1] 8 3 3 3 4 7 2 3 2 5

5 5. Continuous distributions

For a continuous random variable, d*() gives a density, not a point probability. For example, \(P(X=1)=0\) for a continuous variable even though dnorm(1) is positive.

5.3 5.3 Simulation: histogram versus theoretical density

set.seed(140)
normal_sim <- rnorm(2000, mean = 0, sd = 1)

hist(
  normal_sim,
  probability = TRUE,
  breaks = 30,
  col = "grey85",
  border = "white",
  main = "Simulation versus Theoretical Normal Density",
  xlab = "x"
)

curve(
  dnorm(x),
  add = TRUE,
  col = "#D55E00",
  lwd = 3
)

The histogram is empirical; the curve is the theoretical model.

6 6. One reusable classroom dataset: Section A and Section B scores

We now create one labelled data frame and reuse it for the inferential demonstrations.

Suppose two sections take the same examination. Scores are generated only for classroom demonstration.

set.seed(2026)

nA <- 40
nB <- 40

score_A <- round(rnorm(nA, mean = 72, sd = 8), 1)
score_B <- round(rnorm(nB, mean = 76, sd = 9), 1)

# Keep scores within the conventional 0--100 range.
score_A <- pmin(pmax(score_A, 0), 100)
score_B <- pmin(pmax(score_B, 0), 100)

scores <- data.frame(
  ID = sprintf("S%03d", 1:(nA + nB)),
  Section = factor(rep(c("A", "B"), c(nA, nB))),
  Score = c(score_A, score_B)
)

# A categorical outcome will be useful for the chi-square demonstration.
scores$Result <- factor(
  ifelse(scores$Score >= 70, "Pass", "Needs_Improvement")
)

head(scores)
str(scores)
#> 'data.frame':    80 obs. of  4 variables:
#>  $ ID     : chr  "S001" "S002" "S003" "S004" ...
#>  $ Section: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ Score  : num  76.2 63.4 73.1 71.3 66.7 51.9 66.1 63.8 72.9 68.2 ...
#>  $ Result : Factor w/ 2 levels "Needs_Improvement",..: 2 1 2 2 1 1 1 1 2 1 ...
summary(scores)
#>       ID            Section     Score                     Result  
#>  Length:80          A:40    Min.   :51.60   Needs_Improvement:28  
#>  Class :character   B:40    1st Qu.:67.35   Pass             :52  
#>  Mode  :character           Median :73.20                         
#>                             Mean   :73.25                         
#>                             3rd Qu.:77.83                         
#>                             Max.   :95.10
table(scores$Section)
#> 
#>  A  B 
#> 40 40
table(scores$Section, scores$Result)
#>    
#>     Needs_Improvement Pass
#>   A                15   25
#>   B                13   27

The exact sample values are reproducible because the seed is fixed.

6.1 6.1 Visual inspection

par(mfrow = c(1, 2))

boxplot(
  Score ~ Section,
  data = scores,
  col = c("#0072B2", "#D55E00"),
  xlab = "Section",
  ylab = "Score",
  main = "Scores by Section"
)

hist(
  scores$Score[scores$Section == "A"],
  breaks = 10,
  col = adjustcolor("#0072B2", alpha.f = 0.60),
  xlim = range(scores$Score),
  xlab = "Score",
  main = "Section A and B Score Histograms"
)

hist(
  scores$Score[scores$Section == "B"],
  breaks = 10,
  col = adjustcolor("#D55E00", alpha.f = 0.45),
  add = TRUE
)

legend(
  "topright",
  c("Section A", "Section B"),
  fill = c(adjustcolor("#0072B2", .60), adjustcolor("#D55E00", .45)),
  bty = "n"
)

par(mfrow = c(1, 1))

7 7. z-test: known population standard deviation

A z-test is appropriate when the standard deviation used in the standard error is treated as known.

For the one-sample problem

\[ H_0:\mu=\mu_0 \qquad\text{versus}\qquad H_1:\mu\ne\mu_0, \]

the statistic is

\[ Z=\frac{\bar X-\mu_0}{\sigma/\sqrt n}. \]

For demonstration, suppose the historical population standard deviation for Section A scores is known to be \(\sigma=8\), and we test \(H_0:\mu=70\).

A <- scores$Score[scores$Section == "A"]

mu0 <- 70
sigma_known <- 8
alpha <- 0.05

z_stat <- (mean(A) - mu0) / (sigma_known / sqrt(length(A)))
z_pvalue <- 2 * pnorm(-abs(z_stat))

z_critical <- qnorm(1 - alpha / 2)
z_CI <- mean(A) + c(-1, 1) * z_critical *
  sigma_known / sqrt(length(A))

c(
  sample_mean = mean(A),
  z_statistic = z_stat,
  p_value = z_pvalue,
  CI_lower = z_CI[1],
  CI_upper = z_CI[2]
)
#> sample_mean z_statistic     p_value    CI_lower    CI_upper 
#>  71.6725000   1.3222273   0.1860925  69.1933199  74.1516801

7.1 7.1 Using BSDA::z.test()

library(BSDA)

BSDA::z.test(
  A,
  mu = 70,
  sigma.x = 8,
  alternative = "two.sided",
  conf.level = 0.95
)

The output should be read in the same order as any inferential procedure:

estimate → test statistic → p-value → confidence interval → conclusion in context.

8 8. Student’s t distribution and t-tests

When the population standard deviation is unknown and estimated by the sample standard deviation, the standardized mean follows a Student \(t\) distribution under the normal model:

\[ T=\frac{\bar X-\mu_0}{S/\sqrt n}. \]

8.1 8.1 Normal versus t

curve(
  dnorm(x),
  from = -4, to = 4,
  lwd = 3,
  ylab = "Density",
  main = "Normal and Student t Distributions"
)

curve(dt(x, df = 5), add = TRUE, col = "#D55E00", lwd = 3, lty = 2)
curve(dt(x, df = 20), add = TRUE, col = "#0072B2", lwd = 3, lty = 3)

legend(
  "topright",
  c("N(0,1)", "t, df=5", "t, df=20"),
  col = c("black", "#D55E00", "#0072B2"),
  lty = c(1, 2, 3),
  lwd = 3,
  bty = "n"
)

As the degrees of freedom increase, the \(t\) distribution approaches the standard normal distribution.

8.2 8.2 One-sample t-test

t.test(
  A,
  mu = 70,
  alternative = "two.sided",
  conf.level = 0.95
)
#> 
#>  One Sample t-test
#> 
#> data:  A
#> t = 1.3244, df = 39, p-value = 0.1931
#> alternative hypothesis: true mean is not equal to 70
#> 95 percent confidence interval:
#>  69.11814 74.22686
#> sample estimates:
#> mean of x 
#>   71.6725

8.3 8.3 Two-sample t-test using the same data frame

We now test whether the two section means differ:

\[ H_0:\mu_A=\mu_B. \]

The default t.test() uses Welch’s two-sample procedure, which does not require equal population variances.

t.test(
  Score ~ Section,
  data = scores,
  alternative = "two.sided",
  conf.level = 0.95
)
#> 
#>  Welch Two Sample t-test
#> 
#> data:  Score by Section
#> t = -1.7536, df = 77.971, p-value = 0.08342
#> alternative hypothesis: true difference in means between group A and group B is not equal to 0
#> 95 percent confidence interval:
#>  -6.7527786  0.4277786
#> sample estimates:
#> mean in group A mean in group B 
#>         71.6725         74.8350

The confidence interval is for the difference in population means.

9 9. Chi-square distribution and chi-square test

The chi-square distribution appears naturally in inference for variances and in categorical-data tests.

curve(
  dchisq(x, df = 5),
  from = 0,
  to = 20,
  lwd = 3,
  col = "#009E73",
  xlab = expression(chi^2),
  ylab = "Density",
  main = expression(chi^2 * " Distribution, df = 5")
)

qchisq(0.95, df = 5)
#> [1] 11.0705

9.1 9.1 Chi-square test of association

Using the same score data, consider Section and the demonstration variable Result.

score_table <- table(scores$Section, scores$Result)
score_table
#>    
#>     Needs_Improvement Pass
#>   A                15   25
#>   B                13   27
addmargins(score_table)
#>      
#>       Needs_Improvement Pass Sum
#>   A                  15   25  40
#>   B                  13   27  40
#>   Sum                28   52  80

Test

\[ H_0:\text{Section and Result are independent}. \]

chi_result <- chisq.test(score_table)
chi_result
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  score_table
#> X-squared = 0.054945, df = 1, p-value = 0.8147
chi_result$expected
#>    
#>     Needs_Improvement Pass
#>   A                14   26
#>   B                14   26

The expected frequencies are important because the chi-square approximation relies on sufficiently large expected cell counts.

9.2 9.2 Chi-square confidence interval for a normal population variance

For a normal population,

\[ \frac{(n-1)S^2}{\sigma^2}\sim\chi^2_{n-1}. \]

Hence a \(100(1-\alpha)\%\) confidence interval for \(\sigma^2\) is

\[ \left[ \frac{(n-1)S^2}{\chi^2_{1-\alpha/2,n-1}}, \frac{(n-1)S^2}{\chi^2_{\alpha/2,n-1}} \right]. \]

For Section A:

n <- length(A)
s2 <- var(A)
alpha <- 0.05

variance_CI <- c(
  (n - 1) * s2 / qchisq(1 - alpha/2, df = n - 1),
  (n - 1) * s2 / qchisq(alpha/2, df = n - 1)
)

sd_CI <- sqrt(variance_CI)

variance_CI
#> [1]  42.80587 105.17653
sd_CI
#> [1]  6.54262 10.25556

10 10. F distribution and F test

For independent normal samples,

\[ \frac{S_A^2/\sigma_A^2}{S_B^2/\sigma_B^2} \]

leads to an \(F\) distribution. Under

\[ H_0:\sigma_A^2=\sigma_B^2, \]

the observed variance ratio is compared with an \(F\) distribution.

curve(
  df(x, df1 = 10, df2 = 20),
  from = 0,
  to = 4,
  lwd = 3,
  col = "#CC79A7",
  xlab = "F",
  ylab = "Density",
  main = "F Distribution"
)

qf(0.95, df1 = 10, df2 = 20)
#> [1] 2.347878

Using the same Section A/B scores:

var.test(
  Score ~ Section,
  data = scores,
  alternative = "two.sided",
  conf.level = 0.95
)
#> 
#>  F test to compare two variances
#> 
#> data:  Score by Section
#> F = 0.96223, num df = 39, denom df = 39, p-value = 0.9049
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.5089237 1.8193096
#> sample estimates:
#> ratio of variances 
#>          0.9622316

The confidence interval returned by var.test() is for the population variance ratio.

11 11. Type I error, Type II error, and power: a pictorial demonstration

This section should be taught visually.

Consider the upper-tailed z-test

\[ H_0:\mu=\mu_0 \qquad\text{versus}\qquad H_1:\mu>\mu_0, \]

with known \(\sigma\).

The test rejects when

\[ \bar X > c, \qquad c=\mu_0+z_{1-\alpha}\frac{\sigma}{\sqrt n}. \]

Under \(H_0\),

\[ \alpha=P_{H_0}(\bar X>c) \]

is the Type I error probability.

For a particular alternative \(\mu_1>\mu_0\),

\[ \beta(\mu_1)=P_{\mu_1}(\bar X\le c) \]

is the Type II error probability, and

\[ \pi(\mu_1)=1-\beta(\mu_1) \]

is the power.

11.1 11.1 Picture of Type I and Type II errors

mu0 <- 70
mu1 <- 74
sigma <- 8
n <- 40
alpha <- 0.05

se <- sigma / sqrt(n)
critical <- mu0 + qnorm(1 - alpha) * se

xx <- seq(mu0 - 4*se, mu1 + 4*se, length.out = 1200)

f0 <- dnorm(xx, mean = mu0, sd = se)
f1 <- dnorm(xx, mean = mu1, sd = se)

plot(
  xx, f0,
  type = "l",
  lwd = 3,
  col = "#0072B2",
  ylim = c(0, max(f0, f1) * 1.10),
  xlab = expression(bar(X)),
  ylab = "Density",
  main = "Type I Error, Type II Error, and the Rejection Boundary"
)

lines(xx, f1, lwd = 3, col = "#D55E00")

# Type I region: reject H0 when H0 is true.
id_alpha <- xx >= critical
polygon(
  c(critical, xx[id_alpha], max(xx)),
  c(0, f0[id_alpha], 0),
  col = adjustcolor("#0072B2", alpha.f = 0.28),
  border = NA
)

# Type II region: fail to reject H0 when the alternative is true.
id_beta <- xx <= critical
polygon(
  c(min(xx), xx[id_beta], critical),
  c(0, f1[id_beta], 0),
  col = adjustcolor("#D55E00", alpha.f = 0.28),
  border = NA
)

lines(xx, f0, lwd = 3, col = "#0072B2")
lines(xx, f1, lwd = 3, col = "#D55E00")
abline(v = critical, lty = 2, lwd = 2)

alpha_actual <- 1 - pnorm(critical, mean = mu0, sd = se)
beta_actual <- pnorm(critical, mean = mu1, sd = se)

legend(
  "topright",
  legend = c(
    expression(H[0]),
    expression(H[1]),
    paste0("Type I: alpha = ", round(alpha_actual, 3)),
    paste0("Type II: beta = ", round(beta_actual, 3))
  ),
  col = c("#0072B2", "#D55E00", "#0072B2", "#D55E00"),
  lwd = c(3, 3, 8, 8),
  bty = "n"
)

The blue shaded right tail is the probability of rejecting \(H_0\) when \(H_0\) is true. The orange shaded region is the probability of failing to reject \(H_0\) when the displayed alternative is true.

11.2 11.2 Power as the alternative moves away from the null

For the upper-tailed z-test,

\[ \pi(\mu) = P_\mu(\bar X>c) = 1-\Phi\left(\frac{c-\mu}{\sigma/\sqrt n}\right). \]

mu_grid <- seq(65, 85, length.out = 500)

power <- 1 - pnorm(
  critical,
  mean = mu_grid,
  sd = se
)

plot(
  mu_grid,
  power,
  type = "l",
  lwd = 4,
  col = "#009E73",
  ylim = c(0, 1),
  xlab = expression("True mean " * mu),
  ylab = "Power",
  main = "Power Function of a One-Sided z-Test"
)

abline(h = alpha, lty = 2, col = "grey40")
abline(v = mu0, lty = 2, col = "grey40")
abline(h = 1, lty = 3, col = "grey70")

points(mu0, alpha, pch = 19, cex = 1.2)

text(mu0 + 0.5, alpha + 0.06,
     labels = expression(pi(mu[0]) == alpha),
     pos = 4)

text(82, 0.93, "Power approaches 1", pos = 4)

At the null boundary, the rejection probability is \(\alpha\). As the true mean moves farther into the alternative region, the two sampling distributions separate and the power approaches 1.

11.3 11.3 Dynamic classroom animation

The following chunk is intended for interactive execution in RStudio, not for knitting. It moves the alternative mean progressively away from the null and redraws the Type II region and power.

mu0 <- 70
sigma <- 8
n <- 40
alpha <- 0.05

se <- sigma / sqrt(n)
critical <- mu0 + qnorm(1 - alpha) * se

for (mu1 in seq(70, 82, by = 0.25)) {

  xx <- seq(mu0 - 4*se, mu1 + 4*se, length.out = 1000)

  f0 <- dnorm(xx, mu0, se)
  f1 <- dnorm(xx, mu1, se)

  beta <- pnorm(critical, mean = mu1, sd = se)
  power <- 1 - beta

  plot(
    xx, f0,
    type = "l",
    lwd = 3,
    col = "#0072B2",
    ylim = c(0, max(f0, f1) * 1.15),
    xlab = expression(bar(X)),
    ylab = "Density",
    main = paste0(
      "Alternative mean = ", round(mu1, 2),
      "     Power = ", round(power, 3)
    )
  )

  lines(xx, f1, lwd = 3, col = "#D55E00")

  # Fixed Type I rejection region under H0.
  ia <- xx >= critical
  polygon(
    c(critical, xx[ia], max(xx)),
    c(0, f0[ia], 0),
    col = adjustcolor("#0072B2", alpha.f = 0.25),
    border = NA
  )

  # Type II region under the current alternative.
  ib <- xx <= critical
  polygon(
    c(min(xx), xx[ib], critical),
    c(0, f1[ib], 0),
    col = adjustcolor("#D55E00", alpha.f = 0.25),
    border = NA
  )

  lines(xx, f0, lwd = 3, col = "#0072B2")
  lines(xx, f1, lwd = 3, col = "#D55E00")
  abline(v = critical, lty = 2, lwd = 2)

  legend(
    "topright",
    legend = c(
      "Null distribution",
      "Alternative distribution",
      paste0("beta = ", round(beta, 3)),
      paste0("Power = ", round(power, 3))
    ),
    col = c("#0072B2", "#D55E00", "#D55E00", "#009E73"),
    lwd = c(3, 3, 7, 7),
    bty = "n"
  )

  Sys.sleep(0.08)
}

What students should watch: the rejection boundary stays fixed, the null distribution stays fixed, the alternative distribution moves to the right, the Type II region shrinks, and power rises toward 1.

11.4 11.4 Optional: effect of sample size on power

mu_grid <- seq(70, 82, length.out = 400)
sigma <- 8
alpha <- 0.05

plot(
  mu_grid, rep(NA, length(mu_grid)),
  type = "n",
  ylim = c(0, 1),
  xlab = expression("True mean " * mu),
  ylab = "Power",
  main = "Larger Samples Produce Greater Power"
)

n_values <- c(10, 20, 40, 80)
line_types <- 1:4

for (j in seq_along(n_values)) {
  nj <- n_values[j]
  sej <- sigma / sqrt(nj)
  cj <- 70 + qnorm(1 - alpha) * sej

  power_j <- 1 - pnorm(cj, mean = mu_grid, sd = sej)

  lines(
    mu_grid, power_j,
    lwd = 3,
    lty = line_types[j]
  )
}

legend(
  "bottomright",
  legend = paste("n =", n_values),
  lty = line_types,
  lwd = 3,
  bty = "n"
)

12 12. Brownian motion

A standard Brownian motion \(W(t)\) satisfies \(W(0)=0\), has independent increments, and

\[ W(t+\Delta t)-W(t)\sim N(0,\Delta t). \]

Therefore a path can be simulated by generating normal increments and taking cumulative sums.

set.seed(777)

N <- 1000
Tmax <- 1
dt <- Tmax / N

dW <- rnorm(N, mean = 0, sd = sqrt(dt))
W <- c(0, cumsum(dW))
time <- seq(0, Tmax, length.out = N + 1)

plot(
  time, W,
  type = "l",
  lwd = 2,
  col = "#0072B2",
  xlab = "t",
  ylab = "W(t)",
  main = "Simulated Standard Brownian Motion"
)

abline(h = 0, col = "grey70")

13 13. Brownian bridge

From Brownian motion define

\[ B(t)=W(t)-tW(1),\qquad 0\le t\le1. \]

Then

\[ B(0)=B(1)=0. \]

B <- W - time * W[length(W)]

plot(
  time, B,
  type = "l",
  lwd = 2,
  col = "#D55E00",
  xlab = "t",
  ylab = "B(t)",
  main = "Brownian Bridge"
)

abline(h = 0, col = "grey70")

c(B_at_0 = B[1], B_at_1 = B[length(B)])
#> B_at_0 B_at_1 
#>      0      0

The bridge is important because it appears as the limiting process for the centered empirical distribution function under a fully specified continuous null model.

14 14. Kolmogorov statistic and goodness-of-fit

Let \(F_n\) denote the empirical distribution function and \(F_0\) a fully specified continuous null CDF. The one-sample Kolmogorov–Smirnov statistic is

\[ D_n=\sup_x |F_n(x)-F_0(x)|. \]

Under the classical null,

\[ \sqrt{n}\,D_n \Longrightarrow \sup_{0\le t\le1}|B(t)|, \]

where \(B(t)\) is a standard Brownian bridge. The limiting distribution is the Kolmogorov distribution.

14.1 14.1 Visualize empirical and theoretical CDFs

For a clean demonstration, generate a fresh sample from a fully specified \(N(0,1)\) model.

set.seed(909)
x_ks <- rnorm(100, mean = 0, sd = 1)
plot(
  ecdf(x_ks),
  verticals = TRUE,
  do.points = FALSE,
  lwd = 2,
  col = "#0072B2",
  xlab = "x",
  ylab = "CDF",
  main = "Empirical CDF versus N(0,1) CDF"
)

curve(
  pnorm(x, mean = 0, sd = 1),
  add = TRUE,
  col = "#D55E00",
  lwd = 3
)

legend(
  "bottomright",
  c("Empirical CDF", "N(0,1) CDF"),
  col = c("#0072B2", "#D55E00"),
  lwd = c(2, 3),
  bty = "n"
)

14.2 14.2 One-sample KS test

ks.test(
  x_ks,
  "pnorm",
  mean = 0,
  sd = 1
)
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  x_ks
#> D = 0.07725, p-value = 0.5895
#> alternative hypothesis: two-sided

The null hypothesis is that the observations follow the specified \(N(0,1)\) distribution.

Important: Do not estimate mean and sd from the same observations and then interpret the ordinary one-sample KS p-value as the classical normality-test p-value. Parameter estimation changes the null distribution. A Lilliefors-type adjustment or another suitable normality procedure is needed for that different problem.

15 15. The conceptual map of the lecture

The first part of the lecture is unified by R’s distribution grammar:

\[ \boxed{ d \rightarrow \text{density/PMF},\quad p \rightarrow \text{CDF},\quad q \rightarrow \text{quantile},\quad r \rightarrow \text{simulation} } \]

The inferential part is unified by a single Section A/B score dataset:

\[ \boxed{ \text{scores} \rightarrow z \rightarrow t \rightarrow \chi^2 \rightarrow F } \]

The testing ideas are unified visually:

\[ \boxed{ \alpha=\text{Type I error},\qquad \beta=\text{Type II error},\qquad 1-\beta=\text{power} } \]

and the final stochastic-process sequence is

\[ \boxed{ \text{Normal increments} \rightarrow \text{Brownian motion} \rightarrow \text{Brownian bridge} \rightarrow \text{Kolmogorov statistic} \rightarrow \text{KS goodness-of-fit} } \]

16 16. Closing commands students should remember

# Probability distributions
d*()      density / PMF
p*()      CDF / probability
q*()      quantile
r*()      simulation

# Inference
BSDA::z.test()
t.test()
chisq.test()
var.test()
ks.test()

# Empirical distribution
ecdf()

# Simulation
set.seed()
rnorm()
cumsum()

# Graphics
plot()
curve()
hist()
barplot()
lines()
polygon()
legend()

The purpose of the lecture is not to memorize all these commands. It is to see how the same statistical ideas recur across probability, inference, simulation, and stochastic processes.