Simulation & Variability - R Programming

1. Which set of simulation summaries is correct?
   Answer: Bias is $0$, empirical standard deviation is about $1.58$, and RMSE is about $1.41$. When evaluating a simulation, you need to track three distinct quantities — bias, empirical standard deviation, and RMSE — each computed differently, and the key trap here is confusing sample standard deviation (which uses $n-1$) with RMSE (which uses $n$). With theta_hat <- c(8, 9, 10, 11, 12) and theta <- 10, the errors are $(-2, -1, 0, 1, 2)$. Bias is the mean signed error: $\frac{-2-1+0+1+2}{5} = 0$. RMSE is the square root of the mean squared error: $\sqrt{\frac{4+1+0+1+4}{5}} = \sqrt{2} \approx 1.41$. Empirical standard deviation uses R's sd(), which divides by $n-1 = 4$: $\sqrt{\frac{4+1+0+1+4}{4}} = \sqrt{2.5} \approx 1.58$. This makes option C correct — bias $= 0$, SD $\approx 1.58$, RMSE $\approx 1.41$. A is wrong because it assigns RMSE the value $1.58$, which is actually the sample SD. B flips both: it gives SD as $1.41$ (the population-style calculation) and RMSE as $1.58$ (the sample SD value) — a double swap. D is wrong because it reports a nonzero bias; since the estimates are symmetric around $10$, bias is exactly $0$. The key study tip: remember that sd() in R uses $n-1$ (Bessel's correction), while RMSE divides by $n$. These produce different values even from the same data, and exam questions are specifically designed to exploit that confusion.
2. Assuming the simulation design and estimator distribution remain unchanged, what should the analyst expect?
   Answer: The empirical estimator standard deviation should remain near $0.24$, while the Monte Carlo standard error should decrease to about $0.006$. Whenever you see a question mixing "empirical standard deviation of the estimator" with "Monte Carlo standard error," pause and ask yourself: which quantity depends on sample size, and which reflects an intrinsic property of the estimator? The empirical standard deviation (roughly $0.24$) measures the natural spread of the estimator itself — how much individual simulation outputs vary around the true mean. This is a property of the estimator's distribution, not of how many replications you run. Quadrupling replications from $400$ to $1{,}600$ gives you more observations of that same distribution, so the spread stays near $0.24$. The Monte Carlo standard error (MCSE), on the other hand, is the standard error of the mean across replications: $\text{MCSE} = \sigma / \sqrt{n}$. With $n = 1{,}600$, this becomes $0.24/\sqrt{1{,}600} = 0.24/40 = 0.006$. That's exactly what B predicts — the empirical SD stays near $0.24$ while the MCSE drops to $0.006$. A is wrong because it claims both quantities drop, confusing the MCSE reduction with a change in the estimator's own variability. C gets things backwards — it shrinks the empirical SD (which shouldn't change) while keeping the MCSE fixed (which should shrink). D overcalculates the MCSE reduction; going from $400$ to $1{,}600$ replications multiplies $n$ by $4$, so MCSE divides by $\sqrt{4} = 2$, giving $0.006$, not $0.003$. A useful rule of thumb: more replications reduce your uncertainty about the mean (MCSE), but they don't change the estimator's underlying variability (SD). These two quantities answer fundamentally different questions.
3. What is the approximate Monte Carlo standard error of the estimated mean difference?
   Answer: $0.010$, using the paired-difference standard deviation and dividing by $\sqrt{900}$. When comparing two estimators in a simulation study, the key quantity is the variance of the paired difference across replications. Because the two estimators are applied to the same data sets, their results are correlated — and you must account for that correlation when computing the standard deviation of A - B. The variance of a paired difference is $\sigma^2_{A-B} = \sigma^2_A + \sigma^2_B - 2\rho\,\sigma_A\sigma_B$. Plugging in: $0.50^2 + 0.40^2 - 2(0.80)(0.50)(0.40) = 0.25 + 0.16 - 0.32 = 0.09$, so $\sigma_{A-B} = 0.30$. The Monte Carlo standard error of the mean difference is then $0.30 / \sqrt{900} = 0.30/30 = 0.010$. That confirms A is correct. B makes the mistake of assuming independence between estimators, computing $\sqrt{0.25 + 0.16}/\sqrt{900} \approx 0.021$. Ignoring the positive correlation inflates the estimated variance — a common and consequential error when paired structure exists. C uses the correct paired standard deviation ($0.30$) but divides by $\sqrt{100}$ instead of $\sqrt{900}$, as if only 100 replications were used. This inflates the MCSE by a factor of 3, reflecting a miscount of the sample size. D simply reports $\sigma_{A-B} = 0.30$ without dividing by anything, confusing the standard deviation of a single difference with the standard error of the mean difference — a fundamental mix-up between spread and precision. Your study tip: whenever two methods share the same simulated data, always use the paired variance formula. Forgetting the correlation term is one of the most frequent errors in simulation-based comparisons.
4. Which interpretation of the two reported quantiles is most accurate?
   Answer: They describe the central $90\%$ of simulated estimator values, but do not by themselves estimate confidence-interval coverage. When you run a simulation study and call quantile() on your collected estimator values, you're summarizing the empirical distribution of those estimates — nothing more, nothing less. Keep that scope in mind as you evaluate each choice. quantile(theta_hat, c(0.05, 0.95)) returning 1.2 and 2.8 means that 5% of the simulated $\hat{\theta}$ values fell below 1.2 and 5% fell above 2.8. In other words, the central 90% of your estimator's sampling distribution — as approximated by simulation — sits between those two values. That's exactly what D says, and it's the only claim the output directly supports. A is wrong because quantile() says nothing about uncertainty in the mean of theta_hat. A confidence interval for that mean would require a standard error calculation (e.g., using the standard deviation of theta_hat divided by $\sqrt{n}$), not quantiles of the raw estimates. B confuses two different things: the spread of $\hat{\theta}$ values and confidence-interval coverage. Coverage is the proportion of intervals (each built from one data set) that contain the true parameter $\theta = 2$. Computing that requires checking each interval individually — the quantiles of theta_hat don't tell you this. C is a classic trap. The true value $2$ does happen to sit exactly midway between 1.2 and 2.8, but symmetry around the truth in one simulation run doesn't establish unbiasedness. Unbiasedness requires $E[\hat{\theta}] = \theta$, which you'd assess with the mean of theta_hat, not its quantiles. A useful rule of thumb: always ask what object a statistic describes. Quantiles of $\hat{\theta}$ describe $\hat{\theta}$'s distribution — not interval coverage, not means, not bias.
5. Using a normal approximation for Monte Carlo error, which approximate $95\%$ interval describes uncertainty in the estimated coverage probability?
   Answer: $[0.914,\ 0.946]$, obtained from the observed coverage and its binomial Monte Carlo standard error. Whenever you see a simulation study reporting coverage probability, recognize that the estimated coverage is itself a sample proportion — and sample proportions have a well-known binomial standard error that quantifies Monte Carlo uncertainty. Here, the observed coverage is $\hat{p} = 930/1000 = 0.930$. The binomial standard error is $\text{SE} = \sqrt{\hat{p}(1-\hat{p})/n} = \sqrt{(0.930)(0.070)/1000} = \sqrt{0.0000651} \approx 0.00807$. A 95% confidence interval uses $\pm 1.96 \times \text{SE} \approx \pm 0.0158$, giving $[0.930 - 0.016,\ 0.930 + 0.016] \approx [0.914,\ 0.946]$. That confirms C is correct. A is wrong because it centers the interval on the nominal coverage of 0.950 rather than the observed coverage of 0.930. You must use what the simulation actually produced, not the target value. B is wrong because it divides the standard error by an extra $\sqrt{1000}$, effectively computing a standard error of a mean of proportions — a double-shrinkage error. The binomial SE already accounts for the sample size $n = 1000$; there is no reason to divide again. D is wrong because it treats the raw noncoverage count (70) as the standard error, which is dimensionally nonsensical. Standard errors are in probability units (between 0 and 1), not raw counts in the hundreds. Study tip: Anytime a simulation produces a proportion, immediately think "binomial SE = $\sqrt{\hat{p}(1-\hat{p})/n}$." Questions on Monte Carlo error almost always test whether you apply this formula correctly to the observed proportion — not the nominal one.
6. Which report most appropriately summarizes both performance and variability?
   Answer: Report conditional bias $-0.02$ and standard deviation $0.21$ among successful fits, together with a $4\%$ failure rate. When simulation studies encounter convergence failures, the key question is: what population does your summary actually describe? Using na.rm = TRUE silently restricts your analysis to successful replications, so you must be transparent about that restriction and report the failure rate alongside your statistics. Answer B is correct because it does exactly this — it reports the conditional bias ($\hat{\theta} - \theta_0 = 1.98 - 2.00 = -0.02$) and standard deviation ($0.21$) among the 960 successful fits, while explicitly flagging the 4% failure rate. This gives readers a complete, honest picture: the estimator's precision when it works, and how often it doesn't. Answer A is tempting because na.rm = TRUE is syntactically convenient, but it doesn't make failures "irrelevant" — it silently excludes them. Reporting results without acknowledging the 4% failure rate misrepresents the estimator's overall reliability. Answer C sounds like it avoids selection bias by including all 1,000 replications, but substituting failed estimates with 0 introduces severe artificial bias. You're not recovering missing information — you're fabricating it with an arbitrary value that distorts both the mean and standard deviation. Answer D overcorrects in the opposite direction. A 4% failure rate is meaningful, but it doesn't invalidate the information from the 960 successful replications. Discarding valid results entirely is unnecessarily wasteful and tells you nothing about estimator accuracy when convergence succeeds. Study tip: Whenever you see na.rm = TRUE in a simulation context, ask yourself whether the NAs are random noise or systematic failures — if they're systematic, always report the failure rate explicitly alongside your conditional summaries.
7. If the sample size within each simulated data set is increased to $400$, what empirical standard deviation should be expected approximately?
   Answer: $0.15$, because the standard deviation of a sample mean scales approximately as the inverse square root of sample size. Whenever you see a question about the variability of a sample mean, your anchor concept should be the standard error formula: the standard deviation of the sample mean equals $\sigma / \sqrt{n}$. This means variability scales as the inverse square root of sample size — a relationship you must internalize. Starting from a baseline of $0.30$ at $n = 100$, increasing to $n = 400$ multiplies the sample size by $4$. Applying the inverse square root rule: $\sqrt{4} = 2$, so the standard deviation is divided by $2$, giving $0.30 / 2 = 0.15$. That confirms D is correct — the empirical standard deviation you'd expect is approximately $0.15$. A is tempting but wrong. Quadrupling the sample size divides the standard deviation by $\sqrt{4} = 2$, not by $4$. Dividing by $4$ would require sixteenfold the original sample size. This is the classic trap of confusing the square root relationship with a linear one. B is wrong because it conflates two different sources of variability. The number of simulation replications controls how precisely you estimate the empirical standard deviation itself, but the sample size within each data set controls the standard deviation of the sample mean being studied. These are separate quantities. C is wrong in both direction and reasoning. Increasing sample size reduces the standard deviation of the sample mean — it never increases it. C reverses the relationship entirely. As a study tip: always ask yourself "does this scale as $1/n$ or $1/\sqrt{n}$?" For standard deviations of sample means, it's always $1/\sqrt{n}$.
8. Which statement correctly summarizes the simulation?
   Answer: Estimated power is $0.64$, and its Monte Carlo standard error is approximately $0.0096$. When simulating statistical power in R, you're running a test many times under a fixed alternative and tracking how often you reject the null. The proportion of rejections is the estimated power, and because each simulation is an independent Bernoulli trial, you can quantify your simulation's precision using the Monte Carlo standard error (MCSE). Since reject is a logical vector of TRUE/FALSE values, mean(reject) gives the proportion of rejections — that's your estimated power, which here equals 0.64. The MCSE for a proportion estimated from $n$ simulations is $\text{MCSE} = \frac{s}{\sqrt{n}}$, where $s$ is the sample standard deviation. Plugging in the given values: $\frac{0.48}{\sqrt{2500}} = \frac{0.48}{50} = 0.0096$. This confirms B is correct. Choice A confuses the roles of mean and sd, treating sd(reject) = 0.48 as the power estimate and mean(reject) as something else entirely — a straightforward swap of the two statistics. Choice C gets the power right (0.64) but uses the raw standard deviation (0.48) as the MCSE, forgetting to divide by $\sqrt{n}$ — a critical omission that ignores how sample size reduces uncertainty. Choice D computes $1 - 0.64 = 0.36$, which would be the estimated Type II error rate ($\beta$), not power; the MCSE calculation happens to be correct, but it's answering the wrong question. A handy rule: in any simulation study, power = mean(reject) and MCSE = sd(reject) / sqrt(n). Keep these two formulas together and you'll avoid all four traps here.
9. Which interpretation best accounts for these summaries?
   Answer: Rare large losses raise the mean and standard deviation, while the median better represents the loss in a typical replication. Whenever you see a question comparing mean, median, and standard deviation, ask yourself: what happens to each statistic when a distribution has a long tail? That intuition is exactly what this question tests. Here, the data tells a clear story. The median is $1.02$ and the 95th percentile is only $1.30$, meaning 95% of replications have losses below $1.30$. Yet the mean is $2.00$ and the standard deviation is $10.00$ — both dramatically larger. The only way this happens is if a tiny fraction of replications produce catastrophically large losses. Those outliers pull the mean upward (since it sums all values) and inflate the standard deviation (which measures spread from that elevated mean), while leaving the median nearly untouched (since it depends only on rank, not magnitude). This makes C correct: rare extreme losses distort the mean and SD, so the median better describes a typical replication. A is backwards — the mean is more sensitive to extreme values than the median, not less. The median is the resistant statistic here. B confuses a perfectly valid distributional shape with a calculation error. Skewed distributions routinely produce a standard deviation far larger than the median; nothing was computed incorrectly. D conflates distributional skewness with Monte Carlo sampling error. The mean–median gap here reflects the shape of the loss distribution, not the number of replications used. As a study tip: when you see mean $\gg$ median alongside a large SD, immediately think right-skewed distribution with outliers — and remember that the median is your go-to measure of center in that situation.
10. What is the best explanation for the observed lack of variability?
    Answer: Resetting the seed inside every replication generates the same pseudo-random sample each time, so the observed standard deviation is misleading. When working with simulation in R, the placement of set.seed() is everything. A random seed initializes R's pseudo-random number generator at a fixed state — meaning every call to a random function after that seed produces an identical sequence of numbers. If you reset the seed to the same value at the start of each replication, you don't get 500 different samples; you get the same sample repeated 500 times. That's exactly what's happening here: set.seed(9) fires inside the replicate() block, so every iteration draws the identical 25 values from rnorm(25), computes the identical mean, and stores it. When all 500 values in x are the same number, sd(x) is necessarily $0$. Answer A correctly identifies this as the culprit. Answer B is mathematically false. The variance of a sample mean of $n$ standard normal observations is $1/n$, not $0$ — so for $n = 25$, the theoretical standard deviation of the sampling distribution is $1/5 = 0.2$. There is real variability; the simulation just can't reveal it due to the seed issue. Answer C is wrong because replicate() collects individual results into a vector — it performs no averaging whatsoever before returning. Answer D misrepresents what set.seed() does: it does not produce independent samples with matching moments; it produces identical samples, which is the opposite of independence. A practical rule to remember: always place set.seed() outside the replicate() call if you want reproducible but varied simulations. Placing it inside is a subtle bug that silently destroys variability — exactly the kind of trap this exam question is designed to expose.

Start this quiz