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.

enter image description here


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.

  • $\begingroup$ 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). $\endgroup$ – JimB Oct 22 '18 at 3:01
  • $\begingroup$ @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. $\endgroup$ – Jair Taylor Oct 22 '18 at 3:20
  • $\begingroup$ 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. $\endgroup$ – JimB Oct 22 '18 at 4:35
  • $\begingroup$ Are the $X_i$'s independent? $\endgroup$ – StubbornAtom Oct 22 '18 at 4:39
  • $\begingroup$ @StubbornAtom Yes, they are independent. Will edit. $\endgroup$ – 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"}]

True density and beta distribution approximation

  • $\begingroup$ 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. $\endgroup$ – Jair Taylor Oct 22 '18 at 6:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.