beta 分佈隨機變量的 argmax 分佈

讓 $ x_i \sim \text{Beta}(\alpha_i, \beta_i) $ 為了 $ i \in I $ . 讓 $ j = \operatorname*{argmax}_{i \in I} x_i $ （任意斷開關係）。什麼是分佈 $ j $ 按照 $ \alpha $ 和 $ \beta $ ? 除了純採樣之外，有沒有一種有效的計算方法？

當。。。的時候 $ x_i $ 是獨立的 $ 1\le i \le d $ 具有分配功能 $ F_i $ 和密度函數 $ f_i, $ 分別地，機會 $ x_j $ 最大的是（根據分佈函數的定義）

$$ \begin{aligned} \Pr(x_j=\max(x_i,i\in\mathcal I)) &= \Pr(x_1 \le x_j, x_2 \le x_j, \ldots, x_d\le x_j) \ &= E\left[F_1(x_j)F_2(x_j)\cdots F_{j-1}(x_j)(1) F_{j+1}(x_j) \cdots F_d(x_j)\right] \ &= \int_{\mathbb{R}}\left[F_1(x_j)\cdots F_{j-1}(x_j)\ F_{j+1}(x_j) \cdots F_d(x_j)\right]f_j(x_j),\mathrm{d}x_j. \end{aligned} $$

前提是沒有 $ \alpha_i $ 和 $ \beta_i $ 真的很小，這很容易通過數值積分獲得，如下面的R函數beta.argmax所示。（當存在一些微小值的可能性時，將需要更複雜的代碼，因為最高密度的區域可能具有溢出雙精度算術的密度值。實際上，“微小”意味著更接近 $ 0 $ 比 $ 1. $ )

作為它的使用示例，我生成了 $ d=8 $ 值 $ \alpha_i $ 和 $ \beta_i. $

d <- 8
set.seed(17)
alpha <- rexp(d) + 0.1
beta <- rexp(d) + 0.1

然後我計算了概率分佈，並通過模擬 100,000 次迭代對其進行了雙重檢查：

p <- beta.argmax(alpha, beta, stop.on.error=FALSE) # The calculation

x <- matrix(rbeta(d * 1e5, alpha, beta), nrow=d)   # The simulated x_j
p.hat <- tabulate(apply(x, 2, which.max), nbins=d) # Summary of the argmaxes

(signif(rbind(Calculated=p, Simulated=p.hat/sum(p.hat)), 3)) # Comparison
chisq.test(p.hat, p=p)                                       # Formal comparison

輸出是

             [,1]   [,2]    [,3]  [,4]  [,5]   [,6]   [,7]  [,8]
Calculated 0.0247 0.0218 0.00230 0.124 0.451 0.0318 0.0341 0.311
Simulated  0.0245 0.0217 0.00225 0.125 0.451 0.0312 0.0346 0.311

  Chi-squared test for given probabilities

data:  p.hat
X-squared = 2.468, df = 7, p-value = 0.9295

第一個數組中顯示的計算和模擬之間的一致性非常好，隨後的卡方檢驗證實了這一點。

我做了其他測試 $ d $ 像….一樣大 $ 200, $ 保留所有 $ \alpha_i $ 和 $ \beta_i $ 更多 $ 0.5, $ 結果與計算一致。對於更大的值 $ d $ 結果惡化，表明存在數值問題。（我測試到 $ d=500. $ ) 通過提高數值積分中的誤差容限來解決這些問題（以計算時間為代價，達到一分鐘）。

這是代碼。

beta.argmax <- function(alpha, beta, ...) {
 lower <- min(qbeta(1e-9, alpha, beta))
 upper <- max(qbeta(1-1e-9, alpha, beta))
 p <- rep(NA_real_, length(alpha))
 for (i in seq_along(p)) {
   ff <- function(x) dbeta(x, alpha[i], beta[i], log=TRUE)
   f <- Vectorize(function(x) sum(pbeta(x, alpha[-i], beta[-i], log.p=TRUE)))
   h <- function(x) exp(ff(x) + f(x))
   p[i] <- integrate(h, lower, upper, ...)$value
 }
 cat(sum(p), "\n") # Optional check: see how close to 1.000000 the sum is
 p / sum(p)
}

引用自：https://stats.stackexchange.com/questions/548202