Unit 8, lesson 57

A Repeated-Sampling Investigation

How to compare an exact benchmark with a reproducible repeated-sampling simulation, read the empirical distribution and empirical variance of its output, and keep the sample size apart from the replication count.

  • Video 57
  • Simulation Error

In a constructed model, a delay is either 0 or 4 minutes, each with probability one half, independently from one observation to the next. What does the average of four such delays look like across many samples? The binomial distribution gives the exact answer, and a simulation can produce ten thousand averages to set beside it. This lesson carries out both, compares them honestly, and keeps the number of observations in one sample separate from the number of times the whole sample is repeated.

Key terms

Repeated sampling

Definition Fix a sampling procedure, here drawing a random sample X1,…,XnX_1, \dots, X_n from a model and computing a statistic T=g(X1,…,Xn)T = g(X_1, \dots, X_n). Repeated sampling carries out the whole procedure RR times with fresh random inputs, giving values T1,…,TRT_1, \dots, T_R, one per repetition. When the repetitions are independent, T1,…,TRT_1, \dots, T_R are independent and identically distributed, each with the sampling distribution of TT.

In plain words Each repetition draws a whole new sample and records one value of the statistic. Taking one sample and copying its average thousands of times is not repeated sampling, and neither is reusing the same draws in several repetitions.

Example In the delay model each XiX_i equals 00 or 44 minutes with probability 12\tfrac12 each, the observations are independent, the sample size is n=4n = 4 and the statistic is

Xˉ4=X1+X2+X3+X44.\bar{X}_4 = \frac{X_1 + X_2 + X_3 + X_4}{4}.

For this lesson the procedure was actually run R=10000R = 10000 times in Python 3.12 with the generator random.Random(20261014). Each uniform draw from random() below 0.50.5 gave a delay of 44 and every other draw gave 00. Consecutive groups of four among the 40,000 draws formed the 10,000 samples, and each sample gave one recorded mean.

Exact benchmark

Definition An exact benchmark is a quantity derived exactly from the model before the simulation output is examined, such as the distribution of the statistic, its mean and variance, or the probability of an event, and then used as the reference for that output. Agreement with a benchmark, up to simulation error, supports the claim that the program carries out the stated model but does not prove it. It does not check that the model describes a real process.

In plain words Work out what the answer should be, wherever that is possible, before looking at simulated numbers. A mismatch much larger than the standard error then points to a mistake in the program or in the benchmark calculation, and a match is evidence about the program, not about the real-world assumptions behind the model.

Example Let KK be the number of observations equal to 44 in one sample. Then K∼Bin⁡(4,12)K \sim \operatorname{Bin}(4, \tfrac12), the total delay is 4K4K, and so Xˉ4=K\bar{X}_4 = K.

xˉ\bar{x}0011223344
P(Xˉ4=xˉ)P(\bar{X}_4 = \bar{x})116\tfrac{1}{16}416\tfrac{4}{16}616\tfrac{6}{16}416\tfrac{4}{16}116\tfrac{1}{16}

One delay has mean 22 and variance σ2=4\sigma^2 = 4, so E[Xˉ4]=2E[\bar{X}_4] = 2 and Var⁡(Xˉ4)=4/4=1\operatorname{Var}(\bar{X}_4) = 4/4 = 1. The table agrees, since the squared deviations 44, 11, 00, 11 and 44 from 22, weighted by the probabilities, give 4+4+0+4+416=1\tfrac{4 + 4 + 0 + 4 + 4}{16} = 1. Two benchmark event probabilities are P(Xˉ4≥3)=5/16=0.3125P(\bar{X}_4 \ge 3) = 5/16 = 0.3125 and P(1≤Xˉ4≤3)=14/16=7/8P(1 \le \bar{X}_4 \le 3) = 14/16 = 7/8.

Empirical distribution

Definition For recorded values t1,…,tRt_1, \dots, t_R, the empirical distribution assigns to each value vv the relative frequency Nv/RN_v/R, where NvN_v is the number of repetitions rr with tr=vt_r = v. It is a probability distribution, with nonnegative values that add to 11, but it is computed from one finite run. It approximates the sampling distribution and generally differs from it.

In plain words The empirical distribution is the histogram of the output, scaled so that the bars add to one. It varies from one independent run to another, while the exact distribution stays fixed.

Example The recorded counts of the sample means 00 to 44 were 647, 2447, 3782, 2491 and 633, which add to 10,000.

xˉ\bar{x}0011223344
Count647244737822491633
Empirical0.06470.24470.37820.24910.0633
Exact0.06250.250.3750.250.0625

The largest difference is 0.00530.0053, at xˉ=1\bar{x} = 1. That is ordinary for 10,000 repetitions, since a relative frequency estimating 0.250.25 has standard error 0.25⋅0.75/10000≈0.0043\sqrt{0.25 \cdot 0.75/10000} \approx 0.0043.

Empirical variance

Definition The empirical variance of recorded values t1,…,tRt_1, \dots, t_R with average tˉ\bar{t} is the variance of their empirical distribution,

Var⁡^=1R∑r=1R(tr−tˉ)2=1R∑r=1Rtr2−tˉ 2,\begin{aligned} \widehat{\operatorname{Var}} &= \frac{1}{R} \sum_{r=1}^{R} (t_r - \bar{t})^2 \\[0.4em] &= \frac{1}{R} \sum_{r=1}^{R} t_r^2 - \bar{t}^{\,2}, \end{aligned}

with denominator RR. The sample variance with denominator R−1R - 1 is a different statistic, so a report should say which one it uses.

In plain words The empirical variance describes the spread of the recorded output itself. It is computed like any variance, with the relative frequencies in the role of probabilities. For R≥2R \ge 2 the sample variance equals the empirical variance times R/(R−1)R/(R - 1), so for large RR the two are close, and they are equal only when every recorded value is the same.

Example The recorded means have sum 2447+2⋅3782+3⋅2491+4⋅633=200162447 + 2 \cdot 3782 + 3 \cdot 2491 + 4 \cdot 633 = 20016, and their squares have sum 2447+4⋅3782+9⋅2491+16⋅633=501222447 + 4 \cdot 3782 + 9 \cdot 2491 + 16 \cdot 633 = 50122. Dividing by 1000010000,

tˉ=2.0016,Var⁡^=5.0122−2.00162=1.00579744.\begin{aligned} \bar{t} &= 2.0016, \\[0.4em] \widehat{\operatorname{Var}} &= 5.0122 - 2.0016^2 \\[0.4em] &= 1.00579744 . \end{aligned}

The exact values are 22 and 11. The average misses 22 by 0.00160.0016, small next to the standard error 1/10000=0.011/\sqrt{10000} = 0.01 of an average of 10,000 independent sample means. With denominator 99999999 the variance would be about 1.005901.00590.

Replication count

Definition The replication count RR is the number of independent repetitions of the whole experiment. It differs from the sample size nn, the number of observations inside one repetition. The sample size decides which distribution is being studied, for example through Var⁡(Xˉn)=σ2/n\operatorname{Var}(\bar{X}_n) = \sigma^2/n, while RR decides how precisely the simulation estimates that distribution. An estimate of the probability qq of an event about the statistic averages RR independent indicators, one per repetition, so its standard error is q(1−q)/R\sqrt{q(1 - q)/R}.

In plain words More observations per sample make each average less variable. More repetitions estimate the distribution of the averages more precisely. The total number of random draws, nRnR, counts inputs, and when n>1n > 1 putting it in place of RR in a standard error makes the simulation look more precise than it is.

Example The event Xˉ4≥3\bar{X}_4 \ge 3 occurred in 2491+633=31242491 + 633 = 3124 of the 10,000 repetitions, an estimate of 0.31240.3124 against the exact 0.31250.3125. Its standard error is

(5/16)(11/16)10000≈0.0046.\sqrt{\frac{(5/16)(11/16)}{10000}} \approx 0.0046 .

Treating the 40,000 draws as repetitions would give about 0.00230.0023, half the correct value. Changing the sample size changes the target instead. With n=2n = 2 the sample mean takes the values 00, 22 and 44 with probabilities 14\tfrac14, 12\tfrac12 and 14\tfrac14, and its variance is 22 rather than 11.

Check yourself

  1. Use the recorded table to estimate P(Xˉ4=2)P(\bar{X}_4 = 2) and compare it with the exact value. Does the agreement show that delays in a real service follow the 0 or 4 model?
Show the answer

The estimate is 3782/10000=0.37823782/10000 = 0.3782, and the exact probability is 6/16=0.3756/16 = 0.375. The difference 0.00320.0032 is smaller than the standard error

0.375⋅0.62510000≈0.0048,\sqrt{\frac{0.375 \cdot 0.625}{10000}} \approx 0.0048,

so the run is consistent with a program that carries out the model. It does not show that real delays follow the model. The simulation explores consequences of the stated assumptions, and whether real delays take the values 0 and 4 with equal probability, independently, has to be checked with real data. Running the same program with the same seed again would reproduce 3782 exactly, which confirms the computation but is not a new replication.

  1. Change the probability of a 4-minute delay to 14\tfrac14, keeping independence and four observations per sample. Find the exact distribution of Xˉ4\bar{X}_4, its mean and variance, and P(Xˉ4≥3)P(\bar{X}_4 \ge 3). Can the recorded counts serve as a simulation of this new model?
Show the answer

The sample mean still equals the count KK of 4-minute delays, now with K∼Bin⁡(4,14)K \sim \operatorname{Bin}(4, \tfrac14). So

P(Xˉ4=k)=(4k)34−k256P(\bar{X}_4 = k) = \binom{4}{k} \frac{3^{4 - k}}{256}

for k=0,…,4k = 0, \dots, 4, which gives 8181, 108108, 5454, 1212 and 11 over 256256. The mean is 4⋅14=14 \cdot \tfrac14 = 1 and the variance is 4⋅14⋅34=344 \cdot \tfrac14 \cdot \tfrac34 = \tfrac34. The event Xˉ4≥3\bar{X}_4 \ge 3 has probability (12+1)/256=13/256≈0.0508(12 + 1)/256 = 13/256 \approx 0.0508. The recorded counts cannot serve, because they come from the model with probability 12\tfrac12. These answers are exact theoretical values, and a new simulation of the changed model would have to be run and compared with this new benchmark.

  1. A report estimates P(Xˉ4≥3)P(\bar{X}_4 \ge 3) as 0.31240.3124 with standard error 0.3124⋅0.6876/40000≈0.0023\sqrt{0.3124 \cdot 0.6876/40000} \approx 0.0023, arguing that 40,000 random draws were made. What is wrong, and what is the correct plug-in standard error?
Show the answer

The event concerns one sample mean, so there is one indicator per repetition, and the estimate averages R=10000R = 10000 independent indicators. The 40,000 draws are inputs, four for each indicator, and they are not separate trials of the event. The correct plug-in standard error is

0.3124⋅0.687610000≈0.0046,\sqrt{\frac{0.3124 \cdot 0.6876}{10000}} \approx 0.0046,

twice the reported value, because 40000/10000=2\sqrt{40000/10000} = 2. Using the draw count understates the Monte Carlo error.

How this lesson connects