一枚硬幣比另一枚更好的概率是多少？

September 7, 2018

假設我們有兩個有偏見的硬幣C1，並且C2都有不同的轉頭概率。

我們折騰C1 n1時間並獲得H1正面，C2 n2時間並獲得H2正面。我們發現一枚硬幣正面的比率高於另一枚硬幣。

我們可以說一枚硬幣比另一枚更好的概率是多少？（這裡更好意味著轉頭的實際概率更高）。

考慮到兩個硬幣相等的事實，很容易計算進行該觀察的概率。這可以通過Fishers 精確檢驗來完成。鑑於這些觀察

$$ \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} $$

在給定嘗試次數的情況下，硬幣相等時觀察到這些數字的概率 $ 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)!}. $$

但是你要問的是一枚硬幣更好的概率。由於我們爭論對硬幣有多大偏見的看法，我們必須使用**貝葉斯方法來計算結果。請注意，在貝葉斯推理中，術語信念被建模為概率**，並且這兩個術語可以互換使用（s.貝葉斯概率）。我們稱硬幣的概率 $ 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)} $$ 概率密度函數 (pdf) $ f(H_i|p_i,n_i) $ 由二項式概率給出，因為單個嘗試是伯努利實驗： $$ 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) $ .

為了計算 $ f(n_i,H_i) $ 我們使用這樣一個事實，即 pdf 上的積分必須為 1 $ \int_0^1f(p|H_i,n_i)\mathrm dp = 1 $ . 因此，分母將是實現這一目標的不變因素。有一個已知的 pdf 與提名者僅相差一個常數因子，即beta 分佈。因此 $$ 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}. $$

獨立硬幣概率對的 pdf 是 $$ 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). $$

現在我們需要將其整合到以下情況中 $ 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} $$

我無法解析地求解最後一個積分，但可以在插入數字後用計算機進行數值求解。 $ B(\cdot,\cdot) $ 是 beta 函數和 $ B(\cdot;\cdot,\cdot) $ 是不完全 beta 函數。注意 $ \mathbb P(p_1=p_2) = 0 $ 因為 $ p_1 $ 是一個連續變量，從不完全相同 $ p_2 $ .

關於先前的假設 $ f(p_i) $ 並對此發表評論：許多人認為模型的一個很好的替代方法是使用 beta 分佈 $ 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 $ 因此，heads 相當於擁有更多數據。 $ a_i + b_i $ 是如果我們事先將其包括在內，我們將不必進行的投擲次數。

OP表示，這兩個硬幣都存在未知程度的偏差。所以我明白所有的知識都必須從觀察中推斷出來。這就是為什麼我選擇了一個不提供信息的之前，不要將結果偏向於普通硬幣。

所有信息都可以通過以下形式傳達 $ (H_i, n_i) $ 每枚硬幣。缺乏信息性先驗僅意味著需要更多的觀察來確定哪個硬幣更好的概率。

這是R中提供函數的代碼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))
}
你可以畫 $ 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)")
你可能需要install.packages("lattice")先。

可以看出，即使具有統一的先驗和較小的樣本量，當 $ H_1 $ 和 $ H_2 $ 差異足夠大。如果需要，則需要更小的相對差異 $ n_1 $ 和 $ n_2 $ 更大。這是一個情節 $ n_1=100 $ 和 $ n_2=200 $ :

Martijn Weterings建議計算兩者之間差異的後驗概率分佈 $ p_1 $ 和 $ p_2 $ . 這可以通過在集合上集成對的 pdf 來完成 $ 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} $$

同樣，不是我可以分析解決的積分，但 R 代碼將是：
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)
}
我畫了 $ 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 = "-")
你可以計算出概率 $ |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 $ . 但結果往往略有偏差。