2

If we repeatedly roll two fair (6 sided) dice and divide the first result by the second number - what number will we expect to see on average?

I think this can be done by manually enumerating all possible ratios and their corresponding probabilities:

$$E(Z) = \sum_{x=1}^{6} \sum_{y=1}^{6} \frac{x}{y} \cdot \frac{1}{36}$$

$$E(Z) = \frac{1}{36} \left[ \frac{1}{1} + \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \frac{1}{5} + \frac{1}{6} + \frac{2}{1} + \frac{2}{2} + \frac{2}{3} + \frac{2}{4} + \frac{2}{5} + \frac{2}{6} + \frac{3}{1} + \frac{3}{2} + \frac{3}{3} + \frac{3}{4} + \frac{3}{5} + \frac{3}{6} + \frac{4}{1} + \frac{4}{2} + \frac{4}{3} + \frac{4}{4} + \frac{4}{5} + \frac{4}{6} + \frac{5}{1} + \frac{5}{2} + \frac{5}{3} + \frac{5}{4} + \frac{5}{5} + \frac{5}{6} + \frac{6}{1} + \frac{6}{2} + \frac{6}{3} + \frac{6}{4} + \frac{6}{5} + \frac{6}{6} \right] = 1.429$$

I wrote an R function to compute this sum:

compute_E_Z <- function() {
    sum <- 0
    for (x in 1:6) {
        for (y in 1:6) {
            sum <- sum + x/y
        }
    }
    E_Z <- sum * 1/36
    return(E_Z)
}

compute_E_Z() #[1] 1.429167 (exact mean)

Now, I thought I should be able to randomly simulate repeated dice rolls and see that the average ratio is in fact close to 1.429.

In this previous question (How do I find a mean and variance of ratio of two random variables for given mean, variance, and co-variance?), the correct formula for the expected value of two random variables is provided:

Let $f(X,Y)$ be the joint PDF of $X$ and $Y$. Taylor series of $f(X,Y) = \frac{X}{Y}$ around $\mu_X$ and $\mu_Y$ gives $$ \begin{split}f(x,y) &= f(\mu_X,\mu_Y) \\&+ (x-\mu_x)\frac{\partial f(x,y)}{\partial x}\Big\vert_{(x,y)=(\mu_X,\mu_Y)} \\&+ (y-\mu_y)\frac{\partial f(x,y)}{\partial y}\Big\vert_{(x,y)=(\mu_X,\mu_Y)} \\& +\frac{1}{2}(x-\mu_x)^2\frac{\partial^2 f(x,y)}{\partial x^2}\Big\vert_{(x,y)=(\mu_X,\mu_Y)} \\&+\frac{1}{2}(y-\mu_y)^2\frac{\partial^2 f(x,y)}{\partial y^2}\Big\vert_{(x,y)=(\mu_X,\mu_Y)} \\&+ (x-\mu_x)(y-\mu_y)\frac{\partial^2 f(x,y)}{\partial x\partial y}\Big\vert_{(x,y)=(\mu_X,\mu_Y)} \end{split}$$

$$E(\frac{X}{Y}) = E(f(x,y)) \simeq \frac{\mu_X}{\mu_Y} - \frac{C_{XY}}{\mu_Y^2} + \frac{\sigma_Y^2 \mu_X}{\mu_Y^3} + \ldots $$

I then performed this simulation in R while comparing it to the simple average of the dice rolls (incorrect):

library(ggplot2)

results_complex <- numeric(1000) results_simple <- numeric(1000)

E_Z_complex <- function(x, y) { mu_X <- mean(x) mu_Y <- mean(y) C_XY <- cov(x, y) sigma_Y <- sd(y) return(mu_X/mu_Y - C_XY/mu_Y^2 + sigma_Y^2 * mu_X/mu_Y^3) }

E_Z_simple <- function(x, y) { return(mean(x/y)) }

for (i in 1:1000) { # Generate 1000 random values for X and Y X <- sample(1:6, 1000, replace = TRUE) Y <- sample(1:6, 1000, replace = TRUE)

# Calculate E(Z) and store the result
results_complex[i] &lt;- E_Z_complex(X, Y)
results_simple[i] &lt;- E_Z_simple(X, Y)

}

df <- data.frame( iteration = rep(1:1000, 2), E_Z = c(results_complex, results_simple), method = factor(rep(c("Complex", "Simple"), each = 1000)) )

trend_complex <- mean(results_complex) trend_simple <- mean(results_simple)

ggplot(df, aes(x = iteration, y = E_Z, color = method)) + geom_line() + geom_hline(yintercept = 1.42, linetype = "dashed", color = "black", size = 1) + geom_hline(yintercept = trend_complex, linetype = "dashed", color = "red", size = 0.5) + geom_hline(yintercept = trend_simple, linetype = "dashed", color = "blue", size = 0.5) + labs(title = "Simulation of E(Z)", x = "Iteration", y = "E(Z)") + scale_color_discrete(name = "Method", labels = c("Correct Formula", "Incorrect Formula")) + theme_bw() + annotate("text", x = 500, y = 1.42, label = "Actual Mean", hjust = 1)

enter image description here

My Question: According to this simulation, it appears that the correct formula is not coming close to the exact mean. Moreso, the incorrect formula appears to be approaching the exact mean.

Can someone please help me understand why this is happening? Do I need to add more terms to the Taylor Expansion to fix this problem? Or is there some other issue? Perhaps this has to do something for the discrete nature of dice and that Taylor Expansions are not intended for this?

Thanks!

stats_noob
  • 3,112
  • 4
  • 10
  • 36

2 Answers2

5

Assuming the dice are distinguishable (e.g., one is red, the other green), and the ratio is calculated consistently (e.g., the value of the red die $X$ divided by the value of the green die $Y$), then the expectation is $$\operatorname{E}[X/Y] = \frac{1}{36} \cdot \frac{6(6+1)}{2} \cdot H_6^{(1)} = \frac{343}{240},$$ where $H_n^{(m)} = \sum_{k=1}^n \frac{1}{k^m}$ is the generalized $n^{\rm th}$ harmonic number of order $m$.

The variance is $$\begin{align} \operatorname{Var}[X/Y] &= \operatorname{E}[\operatorname{Var}[X/Y \mid Y]] + \operatorname{Var}[\operatorname{E}[X/Y \mid Y]] \\ &= \operatorname{E}\left[\frac{6^2 - 1}{12 Y^2}\right] + \operatorname{Var}\left[\frac{7}{2 Y}\right] \\ &= \frac{35}{12} \operatorname{E}[1/Y^2] + \frac{49}{4} \left(\operatorname{E}[1/Y^2] - \operatorname{E}[1/Y]^2 \right) \\ &= \frac{91}{6} \frac{H_{6}^{(2)}}{6} - \frac{49}{4} \left(\frac{H_6^{(1)}}{6}\right)^2 \\ &= \frac{35819}{20736}. \end{align}$$ There is no need to perform a Taylor series approximation, as these are exact results.

My program was coded in a single Mathematica command for $10^7$ simulations:

{Mean[#], Variance[#]} &[Divide @@ # & /@ RandomInteger[{1, 6}, {10^7, 2}]]

resulted in $$\left\{\frac{857689319}{600000000},\frac{621608824283316239}{359999964000000000}\right\}$$

which after converting to a numeric result is approximately $${1.42948, 1.72669},$$ which agrees with the exact computation. I don't know what your code is doing, but if you are not getting similar values for the mean and variance, you have made a programming error.

heropup
  • 135,869
  • thank you so much for your answer! I wanted to work on this question from an estimation perspective. i.e. if you observe data from this dice rolls, how to correctly estimate the mean – stats_noob Dec 26 '23 at 01:25
3

The simple formula is the correct formula for the statistic. You are looking for the average ratio, and the simple formula is calculating it. I don't see why you think it's "wrong".

The complex formula is an approximation. You're calculating the first three terms of a Taylor series, so if it's not giving the right answer, that implies that the fourth term is significant. The question you link to asks how to get an approximation knowing only mean, variance, and covariance. Implied in this question is that one doesn't have the full joint distribution, and have to make do with these three measures of mean, variance, and covariance. But here you do have the full joint distribution, so there's no need to use the approximation.

BTW, $\mu$ is the population mean. In your code, you're calculating the sample mean, which isn't quite the same. Same for variance and covariance.

Acccumulation
  • 12,210
  • thank you so much for your answer! I wanted to treat this problem as an estimation problem. E.g. suppose there are 100 dice rolls - person1 has only access to the ratios whereas person2 has access to the individual dice rolls. Imagine if the outcome of the second dice was related to the outcome of the first dice (e.g. if the first was 3, there is a 0.5 probability that the second dice will also be 3 and 0.5 probability that any other number). In this case, based on the 100 samples - person1 will not be able to calculate the covariances but person2 can. – stats_noob Dec 26 '23 at 04:29
  • Thus, I thought that the Taylor Approximation approach (person2) would result in better mean and variance estimations than person1, since person1 does not have access to the covariances. – stats_noob Dec 26 '23 at 04:30