Designing a Proof by Simulation

Question

I heard in mathematics there is a way of proving things using a technique called "Proof by Simulation". I have a question about this concept.

Part 1: To introduce this concept, consider the following example.

In statistics, there is a popular misconception about Confidence Intervals. It is commonly and falsely believed that a 95% Confidence Interval created for a parameter estimate using the observed data: the true parameter has a 95% chance (i.e. 0.95 probability) of being contained in this interval. There is wrong on multiple levels: The real value of the parameter is a fixed number so the concept of probability does not apply (i.e. it either is or it is not). But further more, the idea of coverage probability can easily disprove this misconception.

For example, suppose in a coin, the true probability of success is 0.01 - we can flip this coin 100 times, calculate the confidence interval and see if the confidence interval contains 0.01. Repeat this 100 times and see what percent of times the confidence interval on average contains 0.01. Repeat this for all p's between 0 and 1 (by 0.01 increments) and we can quickly "prove by simulation" that this fact about the confidence interval is incorrect. Here is some R code to demonstrate this:

# Initialize 
p_values <- seq(0, 1, 0.01)
coverage <- numeric(length(p_values))
Loop
for (i in 1:length(p_values)) {
    p <- p_values[i]
# Generate samples
samples &lt;- rbinom(100*100, size=1, prob=p)
samples &lt;- matrix(samples, nrow=100, ncol=100)

# Calculate means, variances and 95% confidence intervals
means &lt;- colMeans(samples)
variances &lt;- apply(samples, 2, var)
ci_lower &lt;- means - qnorm(0.975)*sqrt(variances/100)
ci_upper &lt;- means + qnorm(0.975)*sqrt(variances/100)

# Check coverage
coverage[i] &lt;- mean(ci_lower &lt;= p &amp; p &lt;= ci_upper)

}
df <- data.frame(p_values, coverage)

Thus we can see, especially for extreme values of p (i.e. close to 0 and close to 1), the 95% Confidence Interval does not contain the true value 95% of the time. By simulation, we have proved this concept.

Part 2: Now, to illustrate my actual problem.

Suppose there are two dice $X$ and $Y$ (each has six sides with equal probabilities). $X$ is rolled first and then $Y$ is rolled. Here, $X$ and $Y$ can be considered as Random Variables. Let's say we are interested in a function of these Random Variables called $Z$ such that $Z = X/Y$. As we know, ratios of Random Variables are quite complicated to deal with.

Suppose we have 100 pairs of dice rolls, i.e. $(x_1,y_1)$, ... $(x_{100}, y_{100})$.

• Person 1 is given the results of both dice: $(x_1,y_1)$, ... $(x_{100}, y_{100})$
• Person 2 is only given the ratios: $r_1$ = $x_1/y_1$ ..... $r_{100}$ = $x_n/y_n$

In this case, it would be incorrect to simply calculate the mean and variance of Rn directly :

$$E(R) = \frac{1}{n} \sum_{i=1}^{n} R_i$$ $$Var(R) = \frac{1}{n} \sum_{i=1}^{n} (R_i - \bar{R})^2$$

The following link (https://www.stat.cmu.edu/~hseltman/files/ratio.pdf) shows us how to correctly calculate the expected value and variance of $x/y$ :

$$Z1 = E(\frac{X}{Y}) = \frac{\mu_{x}}{\mu_{y}} - \frac{\sigma_{x,y}}{\mu_{x}^2} + \frac{\sigma_{y}^2 \cdot \mu_{x}}{\mu_{y}^3}$$

$$Z2 = Var(\frac{X}{Y}) = \left(\frac{\mu_{x}^2}{\mu_{y}^2}\right) \left(\frac{\sigma_{x}^2}{\mu_{x}^2} - \frac{2 \cdot \sigma_{x,y}}{\mu_{x} \cdot \mu_{y}} + \frac{\sigma_{y}^2}{\mu_{y}^2}\right)$$

It is obvious that mathematically the first way is incorrect, but I was wondering if there is some "proof by simulation" technique that can be used to show that the first way will provide incorrect estimates (e.g. overestimate, underestimate) - the same way as the misconception about Coverage Probability was "disproven" using proof by simulation above.

Part 3: Using R code again, I tried to write a simulation where Person 1 and Person 2 are repeatedly presented with random dice rolls and are tasked with estimating the mean and variance of the ratios:

First I defined the estimation functions for Person 2

calculate_z1 <- function(d1, d2) {
  # Calculate the means
  mean_d1 <- mean(d1)
  mean_d2 <- mean(d2)
Calculate the covariance and variance
cov_d1_d2 <- cov(d1, d2)
  var_d2 <- var(d2)
Calculate z1
z1 <- (mean_d1 / mean_d2) - (cov_d1_d2 / (mean_d1^2)) + (var_d2 * mean_d1 / (mean_d2^3))
return(z1)
}
calculate_z2 <- function(d1, d2) {
Calculate the means
mean_d1 <- mean(d1)
  mean_d2 <- mean(d2)
Calculate the covariance and variances
cov_d1_d2 <- cov(d1, d2)
  var_d1 <- var(d1)
  var_d2 <- var(d2)
Calculate z2
z2 <- (mean_d1^2 / mean_d2^2) * (var_d1 / (mean_d1^2) - (2 * cov_d1_d2) / (mean_d1 * mean_d2) + var_d2 / (mean_d2^2))
return(z2)
}

Then I ran the simulation:

library(dplyr)
library(ggplot2)
set.seed(123)
initialize
x1_values <- numeric()
x2_values <- numeric()
z1_values <- numeric()
z2_values <- numeric()
Loop
for (n in 1:1000) {
    # Simulate dice rolls
    d1 <- sample(1:6, n, replace = TRUE)
    d2 <- sample(1:6, n, replace = TRUE)
# data frame
df &lt;- data.frame(d1 = d1, d2 = d2)

# Add d1/d2 column
df &lt;- df %&gt;%
    mutate(d1_d2 = d1 / d2)

# Calculate x1 and x2 (person 1)
x1 &lt;- mean(df<span class="math-container">$d1_d2, na.rm = TRUE)
x2 &lt;- var(df$</span>d1_d2, na.rm = TRUE)

# Calculate z1 and z2 (person 2)
z1 &lt;- calculate_z1(df<span class="math-container">$d1, df$</span>d2)
z2 &lt;- calculate_z2(df<span class="math-container">$d1, df$</span>d2)

# store results
x1_values &lt;- c(x1_values, x1)
x2_values &lt;- c(x2_values, x2)
z1_values &lt;- c(z1_values, z1)
z2_values &lt;- c(z2_values, z2)

}
create a data frame
plot_df <- data.frame(
    n = 1:1000,
    x1 = x1_values,
    x2 = x2_values,
    z1 = z1_values,
    z2 = z2_values
)
Melt
plot_df <- reshape2::melt(plot_df, id.vars = "n")
legend
levels(plot_df$variable) <- c("person1 mean estimate (incorrect)", "person1 mean variance estimate (incorrect)", "person2 mean estimate (correct)", "person2 variance estimate correct")
Plot
ggplot(plot_df, aes(x = n, y = value, color = variable)) +
    geom_line() +
    labs(x = "n", y = "Value", color = "Variable") +
    theme_minimal()

Based on this graph, as the sample size increases:

• The mean estimates for Person 1 and Person 2 are similar
• The variance estimates for Person 1 and Person 2 are not similar (even as sample size increases)
• The "variance of the variance" estimate for Person 1 fluctuates a lot even when sample size increases - whereas this does not happen for Person 2: This suggests that the estimator for Person 2 is more "stable"?

My Question: I am not sure if this is a convincing simulation to prove my original concept : that simply calculating the variance of individual ratios for random variables is incorrect.

In the Coverage Probability simulation from Part 1, we can clearly see that the "lines on the graph" are sometimes below the 95% line, thus making it very obvious and clear to see that misconception is indeed a misconception.

(Assuming my R code is correct) Somebody could look at my simulation and say "so what!?" - the earlier comment I made on stability is irrelevant. Better stability is not a "litmus test" for showing that one estimator is "more correct" over the another estimator - perhaps some estimators deceivingly produce less variance but are nonetheless incorrect, thus making them falsely appealing but fundamentally incorrect.

In this case, can some mathematical simulation be designed that is as conclusive and clear as the Coverage Probability simulation (from Part 1) which undoubtedly lets you see that Person 1's estimators are incorrect as they clearly deviate from clear reference?

In my simulation, I am still worried that someone could look at this plot and say that my simulation is not indicative and not conclusive that Person 1 is incorrect in this case.

Thanks!

Answer 1

I will only address "Part 1."

Your statistic $X$ is $1/100$ times a $\text{Binomial}(100, p)$ random variable. This is a distribution with a PMF that a computer can calculate. If you can find $u$ such that $P(|X-p| \le u) = 0.95$ according to this scaled binomial distribution that $X$ follows, then $(X-u, X+u)$ is a valid 95% confidence interval for $p$, by definition of $u$.

Your code deviates from this in two ways, which together explain why you do not see exactly 95% coverage in simulations.

(1) Normal approximation. It is common to approximate the distribution of $X$ with $N(p, p(1-p)/n)$. However this approximation is not always good.

The basic approximation generally improves as $n$ increases (at least $20$) and is better when $p$ is not near to $0$ or $1$. Various rules of thumb may be used to decide whether n is large enough, and p is far enough from the extremes of zero or one... —Wikipedia

(2) Using empirical distributions based on sampling from the scaled binomial distribution, instead of the true distribution itself.

Your simulation only simulates $X$ 100 times, which adds an extra layer of error in your results. Note that even for $p$ near $0.5$, you are not getting exactly 95% coverage. Increasing the number of simulations would get closer to the true underlying distribution.

Furthermore, when $p$ is close to $0$ or $1$, the sample variance may be very small (consider an extreme example of $p = 0.01$ and $n=100$, then there is approximately $0.37$ chance of getting no heads, and therefore a sample variance of zero!), leading to a confidence interval that is too narrow. This partially explains why your coverage is so much lower as $p$ is close to $0$ or $1$. More discussion here.