beta 分佈隨機變量的 argmax 分佈

讓 xi∼Beta(αi,βi) 為了 i∈I . 讓 j=argmaxi∈Ixi （任意斷開關係）。什麼是分佈 j 按照 α 和 β ? 除了純採樣之外，有沒有一種有效的計算方法？

當。。。的時候 xi 是獨立的 1≤i≤d 具有分配功能 Fi 和密度函數 fi, 分別地，機會 xj 最大的是（根據分佈函數的定義）

Pr(xj=max(xi,i∈I))=Pr(x1≤xj,x2≤xj,…,xd≤xj) =E[F1(xj)F2(xj)⋯Fj−1(xj)(1)Fj+1(xj)⋯Fd(xj)] =∫R[F1(xj)⋯Fj−1(xj) Fj+1(xj)⋯Fd(xj)]fj(xj),dxj.

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

作為它的使用示例，我生成了 d=8 值 αi 和 β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, 保留所有 αi 和 β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