# How do I find data to show whether a shaved die is really loaded?

I'm making a video about dice, so I went online and bought some loaded dice. The ones I bought are "shaved dice", or "flats", these ones in specific: http://www.amazon.com/gp/product/B008QDJ4RI/ref=oh_details_o04_s00_i00?ie=UTF8&psc=1

I've been doing chi-squared tests of 30 trials at a time with these, and it's really hard to see a bias in any direction. I'm trying to get 95% confidence, but the best I've gotten is 90% (out of 4 tests, 2 of them were at 90% and 2 were very low). I understand it may be a weak effect, but how do I tell with more confidence whether there's actually an effect or not? Do I do more trials? That seems to make the effect even murkier. At what point can I just shrug my shoulders and say, "Well, I guess shaving the die doesn't actually do anything?" Or is 90% good enough?

• Your samples are quite small and the chi-square test is pretty low-power. You need more data than you have. I usually advise at least 200 rolls for a six-sided die and this kind of testing. With a more specific alternative you can lower the sample size. – Glen_b Aug 25 '13 at 10:04
• You can compare the distribution of the chi-square statistics you obtain to the distribution expected from a fair die. (Do not assume this will be a chi-square distribution, because it is not! Simulate this distribution instead.) However, combining all your experimental results and performing a single chi-squared test will be much more powerful. For instance, your method is unlikely to distinguish a 2% increase in the probabilities of one and six even with 100 experiments, whereas a chi-squared test on all 100*30 = 3000 results has around a 99.9% chance of detecting this bias. – whuber Aug 25 '13 at 17:03
• @whuber Well that's frustrating. Everything I've found on the internet for how to test for whether a die is loaded says to use the chi-squared test. For example, this question on math.stackexchange.com where they basically say 30 is the number of rolls you need. Your comment makes sense, but I don't have the knowledge to make a test out of it. – Cassandra Gelvin Aug 25 '13 at 20:16
• I don't see any recommendation on that math thread concerning 30, Cassandra. But what's so frustrating? All I have said is you should simplify your life by not breaking your observations into groups of 30 and doing lots of tests. Just test all the observations at once--and using a chi-squared test is fine for that. (Intuitively you haven't any reasonable hope of detecting small changes in probabilities with only 30 rolls: all observed relative frequencies will be multiples of $1/30 = 3\frac{1}{3}$%, which is poor precision.) – whuber Aug 25 '13 at 20:21
• So am I still doing a chi-squared test? You said not to assume it's a chi-square distribution, but I don't understand what that means. – Cassandra Gelvin Aug 25 '13 at 23:06

What we want to know is how well can a reasonable test (like a chi-squared test) detect a small difference in the chances (of the six outcomes). This is its power: it depends on the size of the difference (a large difference is easy to see) and the number of observations (a large number can detect smaller differences).

The shaving purportedly increases the chance of landing one or six (keeping their chances approximately equal). A natural way to measure the effect of shaving, then, is in terms of the increase in the chance of a one. Let's call this $\varepsilon$. The chances of $(1,2,3,4,5,6)$ are therefore $(1/6+\varepsilon, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6+\varepsilon)$. Clearly $-1/6\le \varepsilon\le 1/3$; here I will examine only cases where $\varepsilon \ge 0$ (shaving favors one and six).

It's not so easy to work out the results of the chi-squared test theoretically in this situation (although it is possible). That's typical of power studies, which usually then resort to Monte-Carlo simulation. Lots of simulations, actually: we need to contemplate various combinations of $\varepsilon$ (which, after all, we don't know at the outset) and sample size $n$. This will give us the information needed to identify the smallest sample size (number of rolls) that will detect an "interesting" difference $\varepsilon$ "reliably."

The individual simulations run and re-run your experiment with $n$ throws of a die, doing this thousands of times. In each experiment a chi-squared statistic is computed. This is a measure of the discrepancy between the distribution of observations and the hypothetical fair distribution ($1/6$ of the total for each face). The proportion of times this statistic is significantly large is the Monte-Carlo estimate of the power of the test (for this particular $n$ for this particular die).

Here, for example, are the outcomes of a set of simulations for $n=30$, setting $\varepsilon$ to a range of values (including $0$, corresponding to a fair die). Each simulation ran $10,000$ times, which is enough to produce robust, reproducible results.

Each simulation is depicted by a histogram of its chi-squared values. The "critical value" for 95% confidence is shown by a vertical red line. All outcomes larger than the critical value are painted in red. In the upper right--the case of no effect--the red values should comprise about 5% of the total, by design, because having 95% confidence means there is a 100-95 = 5% chance of detecting an effect that isn't really there.

As the effect size increases, the chance of detecting it--as indicated by the proportion of red in the histogram--also increases, as expected. It never reaches 100%, though, even for $\varepsilon=1/6=16.7$%. This is an extreme effect: now the one has a $1/6+1/6=1/3$ chance, the six also has a $1/3$ chance, and the other four faces only have a $1/12$ chance apiece. Obviously, tossing the die only $n=30$ times will only be able to detect a gross bias.

Upon repeating the same exercise for a large range of sample sizes $n$, for each possible effect size $\varepsilon$ we can plot the power against $n$. Here are the results:

The same set of effects is represented. We can tell which curve belongs to which effect, because the lower curves (with less power) must correspond to the smaller effects. Thus,

• The bottom (orange) curve is close to $0.05$, especially for large $n$. This is because the tests are run with $100 - 5 =95$% confidence.

• The next lowest curves (light green and teal) show the power for effects of $\varepsilon=1$% and $2.1$%, respectively. Even for sample sizes of $1000$ (that is, $1000$ tosses of the die) we are not likely to detect this much difference in the probabilities. Notice that out of $1000$ tosses with $\varepsilon=2.1$%, we expect there to be about $188$ ones, $188$ sixes, and $156$ each of two through five: that's a pretty big discrepancy for a gambler.

• With just $n=30$ tosses, the only effect we have a reasonable chance of detecting is greater than $8.3$%, where the power is still only about $1/3$.

• (Notice that for smaller sample sizes, less than $30$ or so, the power for $\varepsilon=0$ drops noticeably below the nominal value of $5$%. This is because the chi-squared statistic does not exactly follow a chi-squared distribution for small sample sizes, whence the critical value (which is computed from an assumed chi-squared distribution) is incorrect. If it were to be corrected, the power for the smaller sample sizes would drop a little.)

Evidently, to detect small changes (say of $\varepsilon=1$% or less) many thousands of rolls will be needed. You can work out the power for any $n$ and $\varepsilon$ yourself by modifying and re-running the R code used to produce these figures. Be careful! The study reported here required almost four minutes to run. Test any modifications first on smaller simulations, reducing the number of iterations 1e4 to $1000$ (1e3) or even $100$ at first until you know the code will do exactly what you need.

simulate <- function(eps, n=1) {
# Run a single experiment and return its chi-squared statistic.
d <- sample.int(6, n, prob=rep(1/6,6)+eps*c(1,rep(-1/2,4),1), replace=TRUE)
d <- factor(d, levels=1:6)
chisq.test(tabulate(d))$statistic #$
}
power <- function(eps, n.sample, n.iter) {
# Return the results of n.iter experiments.
replicate(n.iter, simulate(eps, n.sample))
}
power.plot <- function(eps, n.sample, n.iter, alpha=0.05, plot=TRUE) {
# Optionally plot a histogram from power and return its power
# for a test run at 100 - alpha% confidence.
title <- paste("n=", n.sample, " at eps=", round(eps*100, 1), "%", sep="")
d <- power(eps, n.sample, n.iter)
q <- qchisq(1-alpha,5)
if (plot) {
h <- hist(d, plot=FALSE)
hist(d, breaks=h$breaks, main=title, col="#c0c0c0", xlab="Chi-squared", freq=FALSE) h2 <- hist(d[d > q], breaks=h$breaks, plot=FALSE)
h2$density <- h2$density * sum(d>q)/n.iter
plot(h2, add=TRUE, col="#ff404080", freq=FALSE)
abline(v = qchisq(1-alpha,5), col="#ff4040", lwd=2)
}
return (sum(d > q) / length(d))
}
#
# Set up the effects and sample sizes to study.
#
eps <- c(0, 1, 2, 4, 8, 16)/16 * 1/6
n <- c(10, 15, 20, 30, 40, 60, 80, 120, 160, 240, 320, 480, 640, 960)
alpha <- 0.05
#
# Figure 1.
#
par(mfrow=c(2,3))
system.time(p <- sapply(eps, function(eps)
power.plot(eps, n.sample=30, n.iter=1e4, alpha=alpha))
)
#
# Figure 2.
#
system.time(p <- sapply(eps, function(eps)
sapply(n, function(n)
power.plot(eps, n.sample=n, n.iter=1e4, alpha=alpha, plot=FALSE))
))
dimnames(p) <- list(n=n, eps=eps)

colors <- hsv(seq(1/10, 9/10, length.out=length(eps)), .8, 1)
par(mfrow=c(1,1))
plot(c(min(n), max(n)), c(0,1), type="n", log="x",
xlab="Sample size", ylab="Power", main="Power vs. Sample Size")
abline(h=alpha, lwd=2, col="#505050")
tmp <- sapply(1:length(eps),
function(j) lines(n, p[, j], ylim=c(0,1),
col=colors[j], lwd=3, lty=3))

• Wow, when I asked you to make it into an answer you didn't mess around! Thanks so much! Sorry for being whiny. – Cassandra Gelvin Aug 27 '13 at 18:24

One simple approach is to focus on a single face, say six (from the link it looks like this is one of the "flat" faces so should come up less often than 1/6 of the time). Then, if you roll the die $n$ times, you can test the hypothesis that $p = 1/6$ using the test statistic

$$Z = \frac{ \hat p - 1/6}{\sqrt{\frac{ (1/6) \times (5/6) }{n} }}$$

where $\hat p$ is the fraction of rolls that come up six.

A quick calculation tells you that if you roll the die 100 times you will get a p-value < 0.05 (often considered "strong evidence" against the null hypothesis that the die is fair) only if your observed fraction of sixes is more than about 0.24 or less than about 0.1. If you roll the die 1000 times, observed fractions > 0.19 or < 0.145 will yield evidence at the 0.05 level to reject the "fair die" hypothesis.

You can reduce the number of required tosses somewhat (but not much) by counting two shaved faces (e.g., one and six), where the null hypothesis is that these faces will come up 1/6 + 1/6 = 1/3 of the time and the relevant statistic is

$$Z = \frac{ \hat p - 1/3}{\sqrt{\frac{ (1/3) \times (2/3) }{n} }}$$

I suspect that the amount of shaving which can go undetected by the naked eye does not modify roll probabilities substantially, so you might be looking at a lot of throws to test out your dice!

• Okay, so if I have a bunch of dice, what test could I use for fairness that would catch a die like this? Assuming I didn't know this die was in there or how it worked, how would I know which faces to test for? – Cassandra Gelvin Aug 25 '13 at 20:19
• Also, I can tell with the naked eye that it's shaved, if I compare it to another die of the same size. – Cassandra Gelvin Aug 25 '13 at 20:19
• Those are good questions, Cassandra. If you suspect shaving, you might consider pairing the counts for faces and their opposites. If you suspect loading, that is probably a poor idea. All in all, you're probably best off basing your analysis on the finest-grained data you have: the counts of all six individual outcomes. – whuber Aug 26 '13 at 22:58
• Closely related: stats.stackexchange.com/questions/44489/…. – whuber Aug 27 '13 at 15:43