Having a problem with the expectation of the maximum among $n$ independent random variables $ X_1, X_2 \dots X_n$ all ~ the same class of distributions but not necessarily the same mean and other moments (and parameters). We know that, under independence: $$\mathbb{P} \left(\underset{i}{\max }(X_i)\leq a\right)=\prod _k \mathbb{P} (X_k\leq a)$$
I was looking for the result for various distributions of different means (though under independence) and it seems to match the Monte Carlo simulations but the problem is that so far I am failing to get the result with the Gamma distribution without "fudging", while I am getting it under other distributions. I am either too tired and can't see what I am doing wrong, or there is a mistake in the equality above.
Consider the simplified Gamma distribution(1,$\beta$) with CDF $F(a) = 1-e^{-\frac{a}{\beta }}$. So to take the easiest assumption of i.i.d, with $F$ the CDF (cumulative), under the assumption of i.i.d. we can use $F(a)^n$ for the distribution of the maximum. The density: $$f(y)=\frac{\partial F(y)^n}{\partial y}= \frac{n e^{-\frac{y}{\beta }} \left(1-e^{-\frac{y}{\beta }}\right)^{n-1}}{\beta }$$ The expectation of the maximum, $y$ becomes: $$\mathbb{E}(y)= \int_0^{\infty } y\, f(y) \, \mathrm{d}y= \beta \,H_n$$ where $H_n$ is the $n^{th}$ Harmonic number -- corresponding to $\frac{11 \beta }{6}$, $\frac{49 \beta }{20}$ for $n=3$ and $n=6$ respectively. But my Monte Carlo yields results that multiply the above by $\frac{1}{\beta^2}$, namely $\frac{H_n}{\beta}$. What am I doing wrong?
Note 1: I get the same results when using integration by parts, $$\mathbb{E}(y)= \int_0^{\infty } 1-\left(1-e^{-\frac{y}{\beta }}\right)^n \, \mathrm{d}y=\beta \,H_n$$ $\textbf{Solved}$ It turned out to be a bug in Random[GammaDistribution] in Mathematica. I wonder if this post should be deleted.