[Math] Maximum of beta-distributed random variables

Question

Let $X_i \sim \operatorname{Beta}(\alpha_i, \beta_i)$ be independent beta-distributed random variables for $i = 1, \ldots, k$. What can we say about $$X =\max(X_1, \ldots, X_k)?$$ In particular, can we estimate $\alpha$ and $\beta$ so that $X$ is approximately distributed like $\operatorname{Beta}(\alpha, \beta)$? We may assume that $\sum_i \alpha_i + \sum_i \beta_i$ is large if it helps.

We can reduce the question to the case $k =2$, since $$\max(X_1, \ldots, X_k) = \max(\max(\max(X_1, X_2),X_3, \ldots))),$$ although some accuracy might be lost in making successive approximations. Note also that we have $$P(\max(X_1, X_2) \leq z) = P(X_1 \leq z, X_2 \leq z) = P(X_1 \leq z) P(X_2 \leq z),$$ giving the cumulative distribution function for $X$.

Using Sage I was able to take the case $\alpha_1 = 10,\beta_1 = 15,\alpha_2 = 13,\beta_2 = 12$ and approximate the density function of $X$ pretty well with $\alpha = 16.796, \beta = 14.830$. See image.

Context:

This would be useful for the bandit problem or Monte-Carlo tree search. Suppose you are playing $k$ games $Y_i$, and $Y_i$ is either a win, with probability $p_i$, or a loss. Then the game $Y$ which consists of a choice of one of the games $Y_i$ can be modeled by a Bernoulli random variable with parameter $p = \max(p_1, \ldots, p_k)$, since the best strategy is to always choose the game $Y_i$ that has the highest win rate. If we only have limited information about each $Y_i$ (some samples of each, for example), we can put a prior $p_i \sim \operatorname{Beta}(\alpha_i, \beta_i)$ on each parameter $p_i$ and try to infer information about $p$ from this.

Answer 1

One can certainly find a beta distribution with the same mean and variance as $X$ but whether that is a good enough approximation depends on what you need.

If you only want the probability density function of $X$, then that is

$$\sum _{i=1}^n \left(\frac{x^{a_i-1} (1-x)^{b_i-1} \prod _{j \neq i} \frac{B_x(a_j,b_j)}{B(a_j,b_j)}}{B(a_i,b_i)}\right)$$

where $B(a_i,b_i)$ is the beta function and $B_x (a_i,b_i)$ is the incomplete beta function.

I'm not sure there's a nice compact form for the mean and variance but for specific parameters one can calculate the mean and variance which can be matched to a beta distribution. Here's some Mathematica code to do so:

n = 3;
parms = {a[1] -> 1, b[1] -> 6, a[2] -> 4, b[2] -> 7, a[3] -> 4, b[3] -> 5}; 
pdf[x_] := 
  Sum[(x^(a[i] - 1) (1 - x)^(b[i] - 1)/Beta[a[i], b[i]]) Product[Beta[x, a[j], b[j]]/Beta[a[j], b[j]],
   {j, Delete[Range[n], i]}] /. parms, {i, n}];

mean = Integrate[x pdf[x], {x, 0, 1}];
variance = Integrate[x^2 pdf[x], {x, 0, 1}] - mean^2;
sol = N[Solve[{mean == a/(a + b), variance == a b/((a + b)^2 (a + b + 1))}, {a, b}][[1]]]
(* {a -> 6.80319, b -> 6.85957} *)

Plot[{pdf[x], PDF[BetaDistribution[a, b] /. sol, x]}, {x, 0, 1}, 
 PlotLegends -> {"Actual", "Beta approximation"}]