Comprobando si una moneda es justa

Me hizo la siguiente pregunta un amigo. No pude ayudarla, pero espero que alguien me lo explique. No pude encontrar ningún ejemplo similar. Gracias por cualquier ayuda y explicación.

P: Los resultados de 100 experimentos de lanzamiento de monedas se registran como 0 = "Cola" y 1 = "Cabeza". La salida x es una cadena de 0 y 1 de longitud 100. Y el número de veces que obtenemos 1-0-0 en x se calcula y es 20 (ej .: si x = (001001110100), 1-0-0 ocurre 2 veces). ¿Crees que esta es una moneda justa?

Resolviendo el problema por simulación

Mi primer intento sería simular esto en una computadora, que puede lanzar muchas monedas justas muy rápido. A continuación se muestra un ejemplo con un millón de pruebas. El evento 'que el número de veces el patrón' 1-0-0 'ocurre en lanzamientos de monedas es 20 o más' ocurre aproximadamente una vez cada tres mil intentos, por lo que lo que observó no es muy probable (para una feria moneda).$X$$n=100$

Tenga en cuenta que el histrograma es para la simulación y la línea es el cálculo exacto que se explica más adelante.

set.seed(1)

# number of trials
n <- 10^6

# flip coins
q <- matrix(rbinom(100*n, 1, 0.5),n)

# function to compute number of 100 patterns
npattern <- function(x) {
  sum((1-x[-c(99,100)])*(1-x[-c(1,100)])*x[-c(1,2)])
}

# apply function on data 
counts <- sapply(1:n, function(x) npattern(q[x,]))
hist(counts, freq = 0) 

# estimated probability
sum(counts>=20)/10^6
10^6/sum(counts>=20)

---

Resolviendo el problema con un cálculo exacto

Para un enfoque analítico puede usar el hecho de que 'la probabilidad de observar 20 o más secuencias' 1-0-0 'en 100 lanzamientos de monedas es igual a 1 menos la probabilidad de que se necesiten más de 100 lanzamientos para hacer 20 secuencias' . Esto se resuelve en los siguientes pasos:

Tiempo de espera para probabilidad de voltear '1-0-0'

$f_{N,x=1}(n)$

Analicemos las formas de llegar a '1-0-0' como una cadena de Markov. Seguimos los estados descritos por el sufijo de la cadena de volteretas: '1', '1-0' o '1-0-0'. Por ejemplo, si tiene los siguientes ocho lanzamientos 10101100, pasó, en orden, los siguientes ocho estados: '1', '1-0', '1', '1-0', '1', '1', '1-0', '1-0-0' y se necesitaron ocho lanzamientos para llegar a '1-0-0'. Tenga en cuenta que no tiene la misma probabilidad de alcanzar el estado '1-0-0' en cada lanzamiento. Por lo tanto, no puede modelar esto como una distribución binomial . En su lugar, debe seguir un árbol de probabilidades. El estado '1' puede ir a '1' y '1-0', el estado '1-0' puede ir a '1' y '1-0-0', y el estado '1-0-0' es un estado absorbente. Puedes escribirlo como:

           number of flips
           1   2   3   4   5   6   7   8   9   ....   n

'1'        1   1   2   3   5   8  13  21  34   ....   F_n
'1-0'      0   1   1   2   3   5   8  13  21          F_{n-1}
'1-0-0'    0   0   1   2   4   7   12 20  33          sum_{x=1}^{n-2} F_{x}

$n$$n-1$

$$f_{N_c,x=1}(n) = \frac{F_{n-2}}{2^{n-1}}$$

$F_i$$i$

$$f_{N,x=1}(n) = \sum_{k=1}^{n-2} 0.5^{k} f_{N_c,x=1}(1+(n-k)) = 0.5^{n} \sum_{k=1}^{n-2} F_{k}$$

$k$

Esto se puede calcular por convolución.

$$f_{N,x=k}(n) = \sum_{l=1}^{n} f_{N,x=1}(l)f_{N,x=1}(n-l)$$

obtendrá la probabilidad de observar 20 o más patrones '1-0-0' (según la hipótesis de que la moneda es justa)

> # exact computation
> 1-Fx[20]
[1] 0.0003247105
> # estimated from simulation
> sum(counts>=20)/10^6
[1] 0.000337

Aquí está el código R para calcularlo:

# fibonacci numbers
fn <- c(1,1)
for (i in 3:99) {
  fn <- c(fn,fn[i-1]+fn[i-2])
}

# matrix to contain the probabilities
ps <- matrix(rep(0,101*33),33)

# waiting time probabilities to flip one pattern
ps[1,] <- c(0,0,cumsum(fn))/2^(c(1:101))

#convoluting to get the others
for (i in 2:33) {
  for (n in 3:101) {
     for (l in c(1:(n-2))) {
       ps[i,n] = ps[i,n] + ps[1,l]*ps[i-1,n-l]
     }  
  }
}

# cumulative probabilities to get x patterns in n flips
Fx <- 1-rowSums(ps[,1:100])

# probabilities to get x patterns in n flips
fx <- Fx[-1]-Fx[-33]

#plot in the previous histogram
lines(c(1:32)-0.5,fx)

---

Calcular monedas injustas

$x$$n$$p$

Ahora usamos una generalización de los números de Fibonacci:

$$F_{n}(x) = \begin{cases} 1 & \quad \text{if $n=1$} \\ x & \quad \text{if $n=2$} \\ x(F_{n-1}+ F_{n-2}) & \quad \text{if $n>2$} \end{cases}$$

las probabilidades son ahora como:

$$f_{N_c,x=1,p}(n) = (1-p)^{n-1} F_{n-2}((1-p)^{-1}-1)$$

y

$$f_{N,x=1,p}(n) = \sum_{k=1}^{n-2} p(1-p)^{k-1} f_{N_c,x=1,p}(1+n-k) = p(1-p)^{n-1}\sum_{k=1}^{n-2} F_{k}((1-p)^{-1}-1)$$

Cuando trazamos esto, obtienes:

$p=0.5$$p=0.33$

Por lo tanto, si pensabas antes del experimento que la moneda era poco probable injusta, entonces todavía deberías pensar que la moneda es poco probable injusta.

---

$p_{heads} = p_{tails}$

Uno podría mucho más fácilmente probar la probabilidad de una moneda justa contando el número de caras y colas y usar una distribución binomial para modelar estas observaciones y probar si la observación es particular o no.

Sin embargo, puede ser que la moneda esté lanzando, en promedio, un número igual de caras y colas, pero no es justo con respecto a ciertos patrones. Por ejemplo, la moneda podría tener alguna correlación para los lanzamientos sucesivos de monedas (imagino un mecanismo con cavidades dentro del metal de la moneda que se llena con arena que fluirá como un reloj de arena hacia el extremo opuesto del lanzamiento anterior de la moneda, que está cargando la moneda caer más probablemente en el mismo lado que el lado anterior).

$p$

$p=0.45$$p$$p$

# number of trials
set.seed(1)
n <- 10^6

p <- seq(0.3,0.6,0.02)
np <- length(p)
mcounts <- matrix(rep(0,33*np),33)

pb <- txtProgressBar(title = "progress bar", min = 0,
                     max = np, style=3)
for (i in 1:np) {
  # flip first coins
  qfirst <- matrix(rbinom(n, 1, 0.5),n)*2-1
  # flip the changes of the sign of the coin
  qrest <- matrix(rbinom(99*n, 1, p[i]),n)*2-1
  # determining the sign of the coins
  qprod <- t(sapply(1:n, function(x) qfirst[x]*cumprod(qrest[x,])))
  # representing in terms of 1s and 0s
  qcoins <- cbind(qfirst,qprod)*0.5+0.5
  counts <- sapply(1:n, function(x) npattern(qcoins[x,]))

  mcounts[,i] <- sapply(1:33, function(x) sum(counts==x))
  setTxtProgressBar(pb, i)
}
close(pb)

plot(p,colSums(mcounts[c(20:33),]),
     type="l", xlab="p same flip", ylab="counts/million trials", 
     main="observation of 20 or more times '1-0-0' pattern \n for coin with correlated flips")
points(p,colSums(mcounts[c(20:33),]))

---

Usando las matemáticas en estadística

Todo lo anterior está bien, pero no es una respuesta directa a la pregunta.

"¿Crees que esta es una moneda justa?"

Para responder a esa pregunta, uno puede usar las matemáticas anteriores, pero primero debe describir muy bien la situación, los objetivos, la definición de equidad, etc. Sin ningún conocimiento de los antecedentes y las circunstancias, cualquier cálculo será solo un ejercicio matemático y no una respuesta a La pregunta explícita.

Una pregunta abierta es por qué y cómo estamos buscando el patrón '1-0-0'.

• Por ejemplo, tal vez este patrón no era un objetivo, que se decidió antes de hacer la investigación. Tal vez fue algo que 'se destacó' en los datos y fue algo que llamó la atención después del experimento. En ese caso, uno debe considerar que efectivamente está haciendo múltiples comparaciones .
• Otro problema es que la probabilidad calculada anteriormente es un valor p. El significado de un valor p debe considerarse cuidadosamente. Es no la probabilidad de que la moneda es justa. Es, en cambio, la probabilidad de observar un resultado particular si la moneda es justa. Si uno tiene un entorno en el que conoce alguna distribución de la equidad de las monedas, o puede hacer una suposición razonable, entonces puede tener esto en cuenta y usar una expresión bayesiana .
• Lo que es justo, lo que es injusto. Eventualmente, dados suficientes ensayos, uno puede encontrar un poco de injusticia. Pero, ¿es relevante y esa búsqueda no está sesgada? Cuando nos atenemos a un enfoque frecuentista, entonces uno debe describir algo así como un límite por encima del cual consideramos que una moneda es justa (algún tamaño de efecto relevante). Entonces uno podría usar algo similar a la prueba t de dos caras para decidir si la moneda es justa o no (con respecto al patrón '1-0-0').