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))