Maximum of beta-distributed random variables

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.

• How did you obtain the beta distribution parameters for approximating the distribution of the maximum? Using Mathematica and using the method of moments I get $\alpha=16.8681$ and $\beta=14.671$ (although there's not much difference between the two approximations). – JimB Oct 22 '18 at 3:01
• @JimB: I used Sage's find_fit method after computing the density for 100 points. The parameters I get vary a bit depending on where I put the points. – Jair Taylor Oct 22 '18 at 3:20
• If that is a "least squares" approach I don't know why one would want to use such a technique to estimate the parameters of a probability distribution as you don't have a regression situation. – JimB Oct 22 '18 at 4:35
• Are the $X_i$'s independent? – StubbornAtom Oct 22 '18 at 4:39
• @StubbornAtom Yes, they are independent. Will edit. – Jair Taylor Oct 22 '18 at 5:15

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"}]

• Hmm, but this very computationally intensive. I was hoping for something fast, like a simple formula. For the purposes of Monte-Carlo tree search it would need to be computed a very large number of times. – Jair Taylor Oct 22 '18 at 6:05