Computation of distribution parameters for the maximum of two random variables Let $X$ be a beta distributed variable with parameters $a$, $b$. Let $Y$ be a beta distributed variable with parameters $c$, $d$. Let $Z = \max(X, Y)$.
Does anybody know of a fast way to compute the mean and variance of $Z$ given $(a, b, c, d)$?
$X$, $Y$ don't necessarily need to have beta distributions, but they should be similar -- in $[0, 1]$.
What I really want to find is a way to quickly approximate the distribution of the sample max of a bunch of random variables. Let $Z = \max(\{X_1, X_2,\dots , X_n\})$. If I know parameters (not necessarily beta-distribution) for $X_1, X_2,\dots, X_n$, can I quickly compute a set of parameters (same distribution as $X_i$) to form a distribution that closely approximates $Z$?

EDIT:
To clarify, the purpose of this has to do with programming AI for board games. If we are evaluating a position in a chess game and we determine one move leads to a 60% chance of winning while all others have a 20% chance of winning, then we value the position as a 60% win. This is the minimax algorithm... However, what I'm curious about is whether this can be improved upon. In the simple case, we approximate:
$$
\mu_c = \max(\{\mu_1, \mu_2,\dots, \mu_n\})
$$
where $\mu_c$ is the expectation from the current position and $\mu_i$ is the expectation from each of the sub-trees.
It seems like we should be able to do better than this. Can we find a function f that replaces max and closely approximates parameters for the distribution of the sample maximum?
$$
(\mu_c, \sigma_c) = f\left(\{(\mu_1,\sigma_1), (\mu_2,\sigma_2), \dots, (\mu_n, \sigma_n)\}\right)
$$
In this case, we still want $\mu_i$ to be the expectation from each of the nodes. However, $\sigma_i$ does not necessarily need to be variance/standard deviation; just a parameter that represents uncertainty.
Just as:
$$
\max(\{a, b, c\}) = \max\left(\{\max(\{a, b\}), c\}\right)
$$
We should expect:
$$
f\left(\{a, b, c\}\right) \approx f\left(\{f(\{a, b\}), c\}\right)
$$
So $f$ does not necessarily need to take parameters for more than two nodes.
One thing, however, is that $f$ should be relatively fast to compute. If the integrals must be calculated, I think doing double exponential integration might work well, but I'm not sure --  


*

*http://www.johndcook.com/double_exponential_integration.html  

*http://www.codeproject.com/Articles/31550/Fast-Numerical-Integration
 A: I think I can help with this. 
So, when $Z = \max(\{ X_1, X_2, \ldots, X_n\})$ takes on a particular value $Z=z$, this tells you that:


*
* You have one sample $X_i = z$  



*
* Your remaining $n-1$ samples are all $\leq z$


As such, we can say by inspection that
$$
p(Z = z) \propto p(X=z)P(X \leq z)^{n-1}
$$
where I'm using lowercase $p$ to indicate a probability density, and uppercase to indicate an actual probability. $P(X \leq z)$ you can get using the cumulative density function of your distribution. You'll need to find the normalisation term too to ensure the whole thing integrates to $1$;
As for quickly approximating it, I suspect that unless this has a nice closed form, you're going to be looking at numeric methods. But that should be pretty straight forward:


*
* Generate a load of samples using a uniform distribution in the range $[0,1]$



*
* Use importance sampling to weight them by $P(Z=z)$ (remember that when using importance sampling, you do not need a normalised distribution).



*
* Fit your approximating distribution to the data using the importance sampling weights.


and you're done.
Bear in mind that as $n$ grows, most of the probability mass of your distribution will concentrate near $1$. It may be worth using a slightly more complicated sampling procedure than just uniform to ensure you have a good number of samples in this region. 
EDIT: 
It's been pointed out below that I assume identically distributed data in the above. Yeah, that's wrong. I think we're still ok though.
let's take the two variable situation, where we want to find the distribution over $Z = \max(X,Y)$. Just using some standard marginalisation, and assuming independence of $X$ and $Y$, 
$$
p(Z = z) = \int_x \int_y p(Z=z|X=x,Y=y)p(X=x)p(Y=y) dy dx.
$$
(For brevity, I'm not going to keep writing $p(X=x)$, $p(Y=y)$, etc and just use $p(X)$ and $p(Y)$. Hopefully all still clear.)
Given $Z$ is the max of $X$ and $Y$, I'm gonna say
$$
p(Z=z|X,Y) = \delta(\max(X,Y),z)
$$
where the delta function uses its standard definition. So now we have
$$
p(Z = z) = \int_x \int_y \delta(\max(X,Y),z)p(X)p(Y) dy dx.
$$
We can split up the integration ranges into $X>Y$, $Y>X$ and $Y=X$ as follows,
$$
p(Z = z) =  \int_y \int_{x<y}\delta(\max(X,Y),z)p(X)p(Y) dx dy +  \int_x \int_{y<x} \delta(\max(X,Y),z)p(X)p(Y) dy dx + \int_x \int_{y=x} \delta(\max(X,Y),z)p(X)p(Y)dydx.
$$
I think those ranges are correct. They look right on paper.
That makes our max functions nice and easy to evaluate. 
$$
p(Z = z) =  \int_y\int_{x<y} \delta(Y,z)p(X)p(Y) dx dy + \int_x\int_{y<x}  \delta(X,z)p(X)p(Y) dy dx + \int_x \int_{y=x} \delta(X,z)p(X)p(Y)dxdy
$$
which in turn let's us use the sifting property of the delta function to say
$$
p(Z = z) = \int_{x<z} p(X)p(Y=z) dx + \int_{y<z} p(X=z)p(Y) dy + \int_{y=z}p(X=z)p(Y)dy.
$$
Tidying all this up, and using similar notation to before to show probabilites rather than density functions, we get
$$
p(Z = z) = p(Y=z)P(X<z) +  p(X=z)P(Y<z) + 0,
$$
where I'm pretty sure the last term just integrates away to $0$.
Assuming the above is all correct, I'm guessing a similar form applies when dealing with more than 2 variables.
