Lo que queremos saber es cómo una prueba razonable (como una prueba de chi-cuadrado) detectar una pequeña diferencia en las posibilidades (de los seis resultados). Esta es su poder: depende del tamaño de la diferencia (una gran diferencia es fácil de ver) y el número de observaciones (un gran número puede detectar pequeñas diferencias).
El afeitado supuestamente aumenta la probabilidad de aterrizaje one o six (manteniendo sus posibilidades de aproximadamente iguales). Una forma natural para medir el efecto de afeitar, entonces, es en términos del incremento en la probabilidad de que un one. Vamos a llamar a esta $\varepsilon$. Las posibilidades de $(1,2,3,4,5,6)$ por lo tanto $(1/6+\varepsilon, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6-\varepsilon/2, 1/6+\varepsilon)$. Claramente $-1/6\le \varepsilon\le 1/3$; aquí voy a examinar sólo los casos en que $\varepsilon \ge 0$ (afeitado favorece one y six).
No es tan fácil determinar los resultados de la prueba de chi-cuadrado, teóricamente, en esta situación (aunque es posible). Eso es típico de los estudios sobre la energía, que suele recurrir a la simulación Monte-Carlo. Muchas de las simulaciones, realidad: tenemos que contemplar varias combinaciones de $\varepsilon$ (que, después de todo, no sabemos desde el principio) y el tamaño de la muestra $n$. Esto nos dará la información necesaria para identificar el más pequeño tamaño de la muestra (número de rollos) de que va a detectar una "interesante" diferencia " $\varepsilon$ "de forma fiable."
El individuo simulaciones de ejecución y vuelva a ejecutar el experimento con $n$ tiros de un dado, haciendo esto miles de veces. En cada experimento un chi-cuadrado de estadística se calcula. Esta es una medida de la discrepancia entre la distribución de las observaciones y de la hipotética distribución justa ($1/6$ del total para cada cara). La proporción de veces que este dato es significativamente grande es el Monte-Carlo estimación de la potencia de la prueba (para este particular $n$ para este particular morir).
Aquí, por ejemplo, son los resultados de un conjunto de simulaciones para $n=30$, establecimiento $\varepsilon$ a un rango de valores (incluyendo $0$, que corresponde a una feria de morir). Cada simulación corrió $10,000$ veces, lo cual es suficiente para producir robusto, resultados reproducibles.
![Figure 1: histograms]()
Cada simulación se representa mediante un histograma de sus chi-cuadrado de los valores. El "valor crítico" para el 95% de confianza se muestra por una línea vertical de color rojo. Todos los resultados, mayor que el valor crítico están pintadas en rojo. En la parte superior derecha--el caso de no efecto: con el rojo los valores constituyen aproximadamente el 5% del total, por diseño, porque tener una confianza del 95% significa que hay una 100-95 = 5% de probabilidad de detectar un efecto que en realidad no existe.
Como el efecto que aumenta el tamaño, la posibilidad de detectar que--como lo indica la proporción de rojo en el histograma--también aumenta, como se esperaba. Nunca alcanza el 100%, aunque, incluso para $\varepsilon=1/6=16.7$%. Este es un efecto extremo: ahora el one ha $1/6+1/6=1/3$ de probabilidad, el six también tiene un $1/3$ de probabilidad, y las otras cuatro caras sólo tienen un $1/12$ de probabilidad de cada uno. Obviamente, lanzando el dado sólo $n=30$ veces sólo será capaz de detectar un bruto sesgo.
Después de repetir el mismo ejercicio para una amplia gama de tamaños de muestra $n$, para cada posible efecto del tamaño de $\varepsilon$ podemos trazar el poder en contra de la $n$. Aquí están los resultados:
![Figure 2: power curves]()
El mismo conjunto de efectos que se representa. Podemos decir que la curva pertenece a efecto de que, debido a la baja de las curvas (con menos potencia) debe corresponder a los más pequeños efectos. Por lo tanto,
La parte inferior (naranja) de la curva está cerca de a $0.05$, especialmente para las grandes $n$. Esto es debido a que las pruebas se ejecutan con $100 - 5 =95$% de confianza.
La siguiente más baja de las curvas (la luz verde y verde azulado) muestran el poder para los efectos de $\varepsilon=1$$2.1$%%, respectivamente. Incluso para tamaños de muestra de $1000$ (es decir, $1000$ lanzamientos del dado) no es probable para detectar esta diferencia en las probabilidades. Aviso que de $1000$ lanza con $\varepsilon=2.1$%, esperamos que se acerca $188$ ones, $188$ sixesy $156$ cada two través five: eso es una gran discrepancia para un jugador.
Con sólo $n=30$ tiros, el único efecto que tienen una probabilidad razonable de que la detección es mayor que $8.3$%, donde el poder es todavía sólo acerca de $1/3$.
(Observe que para los pequeños tamaños de muestra, a menos de $30$ o así, el poder de $\varepsilon=0$ disminuye considerablemente por debajo del valor nominal de $5$%. Esto es debido a que el chi-cuadrado de estadística no siga exactamente un chi-cuadrado de la distribución de tamaños de muestra pequeños, donde el valor crítico (que se calcula a partir de un supuesto de chi-cuadrado de distribución) es incorrecta. Si se corregirá, el poder de los pequeños tamaños de la muestra iba a caer un poco.)
Evidentemente, para detectar pequeños cambios (digo de $\varepsilon=1$% o menos) muchos miles de rollos serán necesarios. Usted puede trabajar la potencia para cualquier $n$ $\varepsilon$ ti mismo mediante la modificación y volver a ejecutar el R código que se usa para producir estas cifras. Cuidado! El estudio aquí reportado se requiere de casi cuatro minutos a ejecutar. Prueba cualquier modificación de primera en las pequeñas simulaciones, reduciendo el número de iteraciones 1e4 a $1000$ (1e3) o, incluso, $100$ en el primer hasta que usted sepa que el código va a hacer exactamente lo que usted necesita.
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))
