Com que probabilidade uma moeda é melhor que a outra?

8

Digamos que temos duas moedas tendenciosas C1e C2ambas têm probabilidade diferente de virar a cabeça.

Nós jogamos C1 n1tempos e temos H1cabeças, C2 n2tempos e H2cabeças. E descobrimos que a proporção de cabeças para uma moeda é maior que a outra.

Qual é a probabilidade com que podemos dizer que uma moeda é melhor que a outra? (melhor aqui significa maior probabilidade real de virar a cabeça).

Thirupathi Thangavel
fonte
2
Não é um problema de teste de hipótese em que você precisa testar se as estimativas de probabilidade de virar a cara para as moedas são diferentes.
ingênuo
3
Se você deseja uma probabilidade de que uma moeda seja melhor que outra, você pode fazer isso com uma abordagem bayesiana; tudo o que você precisa é uma distribuição prévia da probabilidade. É um cálculo bastante direto.
Glen_b -Reinstala Monica
2
Não há como esses dados fornecerem uma probabilidade de que uma moeda seja melhor que a outra . Você precisa de informações / suposições adicionais sobre uma distribuição anterior de como as moedas são distribuídas. Digamos, numa abordagem bayesiana, que o resultado será muito diferente com base em suas suposições anteriores. Por exemplo, se essas moedas são moedas normais, provavelmente são iguais a priori e você precisará de uma grande quantidade de tentativas para encontrar uma probabilidade decente de que uma é 'melhor' que a outra (porque a probabilidade é alta de que você tem duas moedas semelhantes que testam acidentalmente diferentes) #
Sextus
1
@katosh, é tão exato quanto o anterior. Digamos que você tenha moedas entre p = 0,49 ep = 0,51. Se, nesse caso, você usar a suposição de que as moedas são distribuídas uniformemente entre 0 e 1, as probabilidades que você atribuir serão mais como suposições incorretas atualizadas do que probabilidades atualizadas. Você frequentemente atribui uma alta probabilidade de ser "melhor" à moeda errada. Seria errado considerar o resultado como "a" probabilidade de a moeda ser "melhor" que a outra. Dados inválidos = resultados incorretos. Neste problema, não devemos trabalhar na resolução da matemática, mas na solução do conhecimento.
Sextus Empiricus
1
Dessa forma, é afirmado corretamente. É uma 'crença' e isso não é tão facilmente equiparado a uma 'probabilidade'.
Sextus Empiricus

Respostas:

7

É fácil calcular a probabilidade de fazer essa observação, dado o fato de as duas moedas serem iguais. Isso pode ser feito pelo teste exato de Fishers . Dadas estas observações

coin 1coin 2headsH1H2tailsn1H1n2H2

a probabilidade de observar esses números enquanto as moedas são iguais, dado o número de tentativas , e a quantidade total de cabeças é n1n2H1+H2

p(H1,H2|n1,n2,H1+H2)=(H1+H2)!(n1+n2H1H2)!n1!n2!H1!H2!(n1H1)!(n2H2)!(n1+n2)!.

Mas o que você está pedindo é a probabilidade de que uma moeda seja melhor. Uma vez que discutimos sobre a crença em como as moedas são tendenciosas, temos que usar uma abordagem bayesiana para calcular o resultado. Observe que, na inferência bayesiana, o termo crença é modelado como probabilidade e os dois termos são usados ​​de forma intercambiável (s. Probabilidade bayesiana ). Chamamos a probabilidade de que a moeda jogue cabeças . A distribuição posterior após a observação, para este é dada por Bayes Teorema : A função densidade de probabilidade (pdf)ipipi f ( p i | H i , n i ) = f ( H i | p i , n i ) f ( p i )

f(pi|Hi,ni)=f(Hi|pi,ni)f(pi)f(ni,Hi)
f(Hi|pi,ni)é dada pela probabilidade binomial, já que as tentativas individuais são experimentos de Bernoulli: eu assumir o conhecimento prévio sobre é que pode estar em qualquer lugar entre e com igual probabilidade, portanto, . Portanto, o nominador é .
f(Hi|pi,ni)=(niHi)piHi(1pi)niHi
f(pi)pi01f(pi)=1f(Hi|pi,ni)f(pi)=f(Hi|pi,ni)

Para calcular , usamos o fato de que a integral sobre um pdf deve ser uma . Portanto, o denominador será um fator constante para conseguir exatamente isso. Existe um pdf conhecido que difere do nomeado apenas por um fator constante, que é a distribuição beta . Portanto, f(ni,Hi)01f(p|Hi,ni)dp=1

f(pi|Hi,ni)=1B(Hi+1,niHi+1)piHi(1pi)niHi.

O pdf para o par de probabilidades de moedas independentes é

f(p1,p2|H1,n1,H2,n2)=f(p1|H1,n1)f(p2|H2,n2).

Agora precisamos integrar isso nos casos em que , a fim de descobrir qual a probabilidade da moeda ser melhor do que a moeda : p1>p212

P(p1>p2)=010p1f(p1,p2|H1,n1,H2,n2)dp2dp1=01B(p1;H2+1,n2H2+1)B(H2+1,n2H2+1)f(p1|H1,n1)dp1

Não consigo resolver esta última integral analiticamente, mas é possível resolvê-la numericamente com um computador depois de inserir os números. é a função beta e é a função beta incompleta. Observe que porque é uma variável contínua e nunca é exatamente igual a .B(,)B(;,)P(p1=p2)=0p1p2


Com relação à suposição anterior sobre e observações sobre ela: Uma boa alternativa ao modelo que muitos acreditam é usar uma distribuição beta . Isso levaria a uma probabilidade final Dessa maneira, pode-se modelar um forte viés em relação às moedas regulares por , grande, mas igual . Seria o equivalente a jogar a moeda mais vezes e receber cabeças portanto equivalentes a apenas ter mais dados. é a quantidade de lançamentos que não precisaríamos fazerf(pi)Beta(ai+1,bi+1)

P(p1>p2)=01B(p1;H2+1+a2,n2H2+1+b2)B(H2+1+a2,n2H2+1+b2)f(p1|H1+a1,n1+a1+b1)dp1.
aibiai+biaiai+bi se incluirmos isso antes.

O OP declarou que as duas moedas são tendenciosas em um grau desconhecido. Então eu entendi que todo conhecimento deve ser inferido a partir das observações. É por isso que optei por um desinformativo antes que a dose não viesse o resultado, por exemplo, para moedas comuns.

Toda a informação pode ser transmitida na forma de por moeda. A falta de um prévio informativo significa apenas mais observações para decidir qual moeda é melhor com alta probabilidade.(Hi,ni)


Aqui está o código em R que fornece uma função usando o anterior uniforme : P(n1, H1, n2, H2) =P(p1>p2)f(pi)=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))
}

Você pode desenhar para diferentes resultados experimentais e corrigir , , por exemplo, com este código sniped:P(p1>p2)n1n2n1=n2=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)")

P (p1> p2) para n1 = n2 = 4

Você pode precisar install.packages("lattice")primeiro.

Pode-se ver, que mesmo com a prévia uniforme e um pequeno tamanho da amostra, a probabilidade ou acreditar que uma moeda é melhor pode se tornar bastante sólida, quando e diferem bastante. Uma diferença relativa ainda menor é necessária se e forem ainda maiores. Aqui está uma trama para e :H1H2n1n2n1=100n2=200

insira a descrição da imagem aqui


Martijn Weterings sugeriu calcular a distribuição de probabilidade posterior para a diferença entre e . Isso pode ser feito integrando o pdf do par ao conjunto : p1p2S(d)={(p1,p2)[0,1]2|d=|p1p2|}

f(d|H1,n1,H2,n2)=S(d)f(p1,p2|H1,n1,H2,n2)dγ=01df(p,p+d|H1,n1,H2,n2)dp+d1f(p,pd|H1,n1,H2,n2)dp

Novamente, não é uma integral que eu possa resolver analiticamente, mas o código R seria:

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

Eu tracei para , , e todos os valores de :f(d|n1,H1,n2,H2)n1=4H1=3n2=4H2

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

insira a descrição da imagem aqui

Você pode calcular a probabilidade deestar acima de um valor por . Lembre-se de que a dupla aplicação da integral numérica vem com algum erro numérico. Por exemplo, deve sempre ser igual a pois sempre assume um valor entre e . Mas o resultado geralmente se desvia um pouco.|p1p2|dintegrate(fd, d, 1, n1, H1, n2, H2)integrate(fd, 0, 1, n1, H1, n2, H2)1d01

katosh
fonte
Eu não sei o valor real de p1
Thirupathi Thangavel
1
Sinto muito pela minha notação ruim 😅, mas você não precisa conectar um valor fixo para . O (agora alterado) na integral é a variável na qual você integra. Assim como você pode integrar sem ter um valor fixo para . p1p101x2dxx
Katrash9
1
O teste exato de Fisher é mais especificamente sobre a hipótese de que as moedas têm a mesma probabilidade e os totais marginais são fixos . Este não é o caso neste problema de moeda. Se você fizer o teste novamente, poderá observar outra quantidade total de cabeças.
Sextus Empiricus 09/09
@MartijnWeterings no meu caso, a probabilidade de virar a cabeça para uma moeda é sempre fixa. Isso não é suficiente?
Thirupathi Thangavel
1
@ThirupathiThangavel o problema com o teste de Fisher é sobre os totais marginais não fixos. O modelo de teste exato assume que a probabilidade de cabeças é a mesma e fixa, mas também os marginais são fixados antes do experimento. Essa segunda parte não é o caso. Isso fornece uma probabilidade condicional diferente para valores extremos. No geral, o teste de Fisher será conservador. A probabilidade do resultado dada a hipótese TRUE (ou seja, probabilidade fixa e similar para cabeças, mas não necessariamente totais marginais fixos) é menor que o calculado (você obtém valores p mais altos).
Sextus Empiricus
2

Fiz uma simulação numérica com R, provavelmente você está procurando uma resposta analítica, mas achei que isso poderia ser interessante para compartilhar.

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

Neste primeiro código, eu simplesmente simulo dois sorteios. Aqui você pode ver, é claro theta1 > theta2, que a frequência H1será maior que H2. Observe o diferente N1, N2os tamanhos.

Vamos ver o que podemos fazer com diferentes thetas. Observe que o código não é o ideal. Em 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

Então 0,93 é a frequência dos tempos (de um 100) em que a primeira moeda tinha mais cabeças. Isso parece bom, olhando theta1e theta2usado.

Vamos ver com dois vetores 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")

Lembre-se de que res_vsão as frequências em que H1 > H2, de 100 simulações.

insira a descrição da imagem aqui

Assim, à medida que theta1aumenta, a probabilidade de H1ser maior aumenta, é claro.

Já fiz algumas outras simulações e parece que os tamanhos N1, N2são menos importantes.

Se você estiver familiarizado, Ruse esse código para esclarecer um pouco o problema. Estou ciente de que essa não é uma análise completa e pode ser melhorada.

RLave
fonte
2
Interessante como res_vmuda continuamente quando os thetas se encontram. Eu entendi a pergunta enquanto perguntava sobre o viés intrínseco das moedas depois de fazer apenas uma única observação. Você parece responder que observações alguém faria depois de conhecer o viés.
katosh