¿Con qué probabilidad una moneda es mejor que la otra?

Digamos que tenemos dos monedas sesgadas C1y C2ambas tienen diferentes probabilidades de girar la cabeza.

Lanzamos C1 n1tiempos y obtenemos H1cabezas, C2 n2tiempos y obtenemos H2cabezas. Y encontramos que la proporción de caras para una moneda es más alta que la otra.

¿Cuál es la probabilidad con la que podemos decir que una moneda es mejor que la otra? (mejor aquí significa mayor probabilidad real de girar la cabeza).

Es fácil calcular la probabilidad de hacer esa observación, dado el hecho de que las dos monedas son iguales. Esto puede hacerse mediante la prueba exacta de Fishers . Dadas estas observaciones

$$\begin{array} {r|c|c} &\text{coin }1 &\text{coin }2 \\ \hline \text{heads} &H_1 &H_2\\ \hline \text{tails} &n_1-H_1 &n_2-H_2\\\end{array}$$

La probabilidad de observar estos números mientras las monedas son iguales dada la cantidad de intentos , y la cantidad total de es $n_1$$n_2$$H_1+H_2$

$$p(H_1, H_2|n_1, n_2, H_1+H_2) = \frac{(H_1+H_2)!(n_1+n_2-H_1-H_2)!n_1!n_2!}{H_1!H_2!(n_1-H_1)!(n_2-H_2)!(n_1+n_2)!}.$$

Pero lo que está pidiendo es la probabilidad de que una moneda sea mejor. Dado que discutimos acerca de una creencia sobre cuán sesgadas están las monedas, tenemos que usar un enfoque bayesiano para calcular el resultado. Tenga en cuenta que en la inferencia bayesiana el término creencia se modela como probabilidad y los dos términos se usan indistintamente (s. Probabilidad bayesiana ). Llamamos a la probabilidad de que la moneda arroje . La distribución posterior después de la observación, para este está dada por teorema de Bayes : El función de densidad de probabilidad (pdf)$i$$p_i$$p_i$

$$f(p_i|H_i,n_i)= \frac{f(H_i|p_i,n_i)f(p_i)}{f(n_i,H_i)}$$ $f(H_i|p_i,n_i)$está dada por la probabilidad Binomial, ya que los intentos individuales son experimentos de Bernoulli: I Asumo El conocimiento previo sobre es que podría estar entre y con igual probabilidad, por lo tanto . Entonces el nominador es . $$f(H_i|p_i,n_i) = \binom{n_i}{H_i}p_i^{H_i}(1-p_i)^{n_i-H_i}$$ $f(p_i)$$p_i$$0$$1$$f(p_i) = 1$$f(H_i|p_i,n_i)f(p_i)= f(H_i|p_i,n_i)$

Para calcular usamos el hecho de que la integral sobre un pdf tiene que ser uno . Entonces el denominador será un factor constante para lograr eso. Existe un pdf conocido que difiere del nominador solo por un factor constante, que es la distribución beta . Por lo tanto, $f(n_i,H_i)$$\int_0^1f(p|H_i,n_i)\mathrm dp = 1$

$$f(p_i|H_i,n_i) = \frac{1}{B(H_i+1, n_i-H_i+1)}p_i^{H_i}(1-p_i)^{n_i-H_i}.$$

El pdf para el par de probabilidades de monedas independientes es

$$f(p_1,p_2|H_1,n_1,H_2,n_2) = f(p_1|H_1,n_1)f(p_2|H_2,n_2).$$

Ahora necesitamos integrar esto en los casos en que para descubrir cómo es probable que la moneda sea ​​mejor que la moneda : $p_1>p_2$$1$$2$

$$\begin{align} \mathbb P(p_1>p_2) &= \int_0^1 \int_0^{p‘_1} f(p‘_1,p‘_2|H_1,n_1,H_2,n_2)\mathrm dp‘_2 \mathrm dp‘_1\\ &=\int_0^1 \frac{B(p‘_1;H_2+1,n_2-H_2+1)}{B(H_2+1,n_2-H_2+1)} f(p‘_1|H_1,n_1)\mathrm dp‘_1 \end{align}$$

No puedo resolver esta última integral analíticamente, pero uno puede resolverla numéricamente con una computadora después de conectar los números. es la función beta y es la función beta incompleta. Tenga en cuenta que porque es una variable continua y nunca es exactamente lo mismo que .$B(\cdot,\cdot)$$B(\cdot;\cdot,\cdot)$$\mathbb P(p_1=p_2) = 0$$p_1$$p_2$

---

Con respecto a la suposición anterior sobre y sus comentarios: Una buena alternativa para modelar muchos cree es usar una distribución beta . Esto llevaría a una probabilidad final De esa manera, uno podría modelar un fuerte sesgo hacia las monedas regulares por grandes pero iguales , . Sería equivalente a lanzar la moneda veces adicionales y recibir , por lo tanto, equivale a tener más datos. es la cantidad de lanzamientos que no tendríamos que hacer$f(p_i)$$Beta(a_i+1,b_i+1)$

$$\mathbb P(p_1>p_2) =\int_0^1 \frac{B(p‘_1;H_2+1+a_2,n_2-H_2+1+b_2)}{B(H_2+1+a_2,n_2-H_2+1+b_2)} f(p‘_1|H_1+a_1,n_1+a_1+b_1)\mathrm dp‘_1.$$ $a_i$$b_i$$a_i+b_i$$a_i$$a_i + b_i$ si incluimos esto antes.

El OP declaró que las dos monedas están sesgadas en un grado desconocido. Entonces entendí que todo el conocimiento tiene que inferirse de las observaciones. Es por eso que opté por un poco informativo antes de que la dosis no sesgue el resultado, por ejemplo, hacia las monedas normales.

Toda la información se puede transmitir en forma de por moneda. La falta de un previo informativo solo significa que se necesitan más observaciones para decidir qué moneda es mejor con alta probabilidad.$(H_i, n_i)$

---

Aquí está el código en R que proporciona una función usando el uniforme anterior : P(n1, H1, n2, H2) $=\mathbb P(p_1>p_2)$$f(p_i)=1$

mp <- function(p1, n1, H1, n2, H2) {
    f1 <- pbeta(p1, H2 + 1, n2 - H2 + 1)
    f2 <- dbeta(p1, H1 + 1, n1 - H1 + 1)
    return(f1 * f2)
}

P <- function(n1, H1, n2, H2) {
    return(integrate(mp, 0, 1, n1, H1, n2, H2))
}

Puede dibujar para diferentes resultados experimentales y , , ejemplo, con este código franqueado:$P(p_1>p_2)$$n_1$$n_2$$n_1=n_2=4$

library(lattice)
n1 <- 4
n2 <- 4
heads <- expand.grid(H1 = 0:n1, H2 = 0:n2)
heads$P <- apply(heads, 1, function(H) P(n1, H[1], n2, H[2])$value)
levelplot(P ~ H1 + H2, heads, main = "P(p1 > p2)")

Es posible que necesite install.packages("lattice")primero.

Se puede ver que, incluso con el uniforme anterior y un tamaño de muestra pequeño, la probabilidad o creer que una moneda es mejor puede volverse bastante sólida, cuando y difieren lo suficiente. Se necesita una diferencia relativa aún menor si y son aún mayores. Aquí hay una gráfica para y :$H_1$$H_2$$n_1$$n_2$$n_1=100$$n_2=200$

---

Martijn Weterings sugirió calcular la distribución de probabilidad posterior para la diferencia entre y . Esto se puede hacer integrando el pdf del par sobre el conjunto : $p_1$$p_2$$S(d)=\{(p_1,p_2)\in[0,1]^2|d=|p_1-p_2|\}$

$$\begin{align} f(d|H_1,n_1,H_2,n_2) &= \int_{S(d)}f(p_1,p_2|H_1,n_1,H_2,n_2) \mathrm d\gamma\\ &= \int_0^{1-d} f(p,p+d|H_1,n_1,H_2,n_2) \mathrm dp + \int_d^1 f(p,p-d|H_1,n_1,H_2,n_2) \mathrm dp\\ \end{align}$$

De nuevo, no es una integral que pueda resolver analíticamente, pero el código R sería:

d1 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p + d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

d2 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p - d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

fd <- function(d, n1, H1, n2, H2) {
    if (d==1) return(0)
    s1 <- integrate(d1, 0, 1-d, d, n1, H1, n2, H2)
    s2 <- integrate(d2, d, 1, d, n1, H1, n2, H2)
    return(s1$value + s2$value)
}

He trazado para , , y todos los valores de :$f(d|n_1,H_1,n_2,H_2)$$n_1=4$$H_1=3$$n_2=4$$H_2$

n1 <- 4
n2 <- 4
H1 <- 3
d <- seq(0, 1, length = 500)

get_f <- function(H2) sapply(d, fd, n1, H1, n2, H2)
dat <- sapply(0:n2, get_f)

matplot(d, dat, type = "l", ylab = "Density",
        main = "f(d | 4, 3, 4, H2)")
legend("topright", legend = paste("H2 =", 0:n2),
       col = 1:(n2 + 1), pch = "-")

Puede calcular la probabilidad deestar por encima de un valor por . Tenga en cuenta que la doble aplicación de la integral numérica viene con algún error numérico. Por ejemplo , siempre debe ser igual a ya que siempre toma un valor entre y . Pero el resultado a menudo se desvía ligeramente.$|p_1-p_2|$$d$integrate(fd, d, 1, n1, H1, n2, H2)integrate(fd, 0, 1, n1, H1, n2, H2)$1$$d$$0$$1$

He hecho una simulación numérica con R, probablemente estás buscando una respuesta analítica, pero pensé que esto podría ser interesante para compartir.

set.seed(123)
# coin 1
N1 = 20
theta1 = 0.7

toss1 <- rbinom(n = N1, size = 1, prob = theta1)

# coin 2
N2 = 25
theta2 = 0.5

toss2 <- rbinom(n = N2, size = 1, prob = theta2)

# frequency
sum(toss1)/N1 # [1] 0.65
sum(toss2)/N2 # [1] 0.52

En este primer código, simplemente simulo dos lanzamientos de monedas. Aquí puede ver, por supuesto, que theta1 > theta2, entonces, por supuesto, la frecuencia de H1será mayor que H2. Tenga en cuenta el diferente N1, N2tamaños.

Veamos qué podemos hacer con diferentes thetas. Tenga en cuenta que el código no es óptimo. En absoluto.

simulation <- function(N1, N2, theta1, theta2, nsim = 100) {
  count1 <- count2 <- 0

  for (i in 1:nsim) {
    toss1 <- rbinom(n = N1, size = 1, prob = theta1)
    toss2 <- rbinom(n = N2, size = 1, prob = theta2)

    if (sum(toss1)/N1 > sum(toss2)/N2) {count1 = count1 + 1} 
    #if (sum(toss1)/N1 < sum(toss2)/N2) {count2 = count2 + 1} 
  }

  count1/nsim

}
set.seed(123)
simulation(20, 25, 0.7, 0.5, 100)
#[1] 0.93

Entonces, 0.93 es la frecuencia de veces (de un 100) que la primera moneda tenía más caras. Esto parece estar bien, mirando theta1y theta2usado.

Veamos con dos vectores de thetas.

theta1_v <- seq(from = 0.1, to = 0.9, by = 0.1)
theta2_v <- seq(from = 0.9, to = 0.1, by = -0.1)

res_v <- c()
for (i in 1:length(theta1_v)) {

  res <- simulation(1000, 1500, theta1_v[i], theta2_v[i], 100)
  res_v[i] <- res

}

plot(theta1_v, res_v, type = "l")

Recuerde que res_vson las frecuencias donde H1 > H2, de 100 simulaciones.

Entonces, a medida que theta1aumenta, la probabilidad de H1ser mayor aumenta, por supuesto.

He hecho algunas otras simulaciones y parece que los tamaños N1, N2son menos importantes.

Si está familiarizado R, puede usar este código para arrojar algo de luz sobre el problema. Soy consciente de que este no es un análisis completo, y se puede mejorar.