Expectation of maximum of arithmetic means of i.i.d. exponential random variables Given the sequence $(X_n), n=1,2,... $, of iid exponential random variables with parameter $1$, define:
$$ M_n := \max \left\{ X_1, \frac{X_1+X_2}{2}, ...,\frac{X_1+\dots+X_n}{n} \right\} $$
I want to calculate $\mathbb{E}(M_n)$. Running a simulation leads me to believe that
$$ \mathbb{E}(M_n)=1+\frac{1}{2^2}+\cdots+\frac{1}{n^2} = H_n^{(2)}.$$
Is this correct? If yes, how would one go proving it? I tried using induction and the fact that $M_{n+1}=\max \{M_n, \frac{1}{n}(X_1+\cdots+X_{n+1}) \}$ along with the equality $E(X_1|X_1+\cdots+X_{n+1})=\frac{1}{n}(X_1+\cdots+X_{n+1})$ but didn't manage to accomplish anything.
 A: For any $x>0$ and $n>1$, the following relation holds (with $\mathbb P(M_1\leqslant x)=1-e^{-x}$): 

$$
\mathbb P(M_n \leqslant x)=\mathbb P(M_{n-1} \leqslant x) - e^{-nx}\frac{x^{n-1}n^{n-2}}{(n-1)!}\tag{1}
$$ 

Consequently, $\mathbb P(M_n \leqslant x) = 1 - \sum\limits_{r=1}^{n} e^{-rx} \frac{x^{r-1}r^{r-2}}{(r-1)!}$. Therefore, 
$$\mathbb E[M_n]=\int\limits_{0}^{\infty} \mathbb P(M_n>x) \mathrm dx = \sum\limits_{r=1}^{n}\int\limits_{0}^{\infty}e^{-rx}\frac{x^{r-1}r^{r-2}}{(r-1)!}\mathrm dx = \sum\limits_{r=1}^{n}\frac{1}{r^2}\,.$$

 Proof of  $(1)$:
$$
\mathbb P(M_{n-1} \leqslant x) - \mathbb P(M_n \leqslant x) = e^{-nx}\int\limits_{0}^{x}\int\limits_{0}^{2x-x_1}\ldots\int\limits_{0}^{(n-1)x-\sum_{i=1}^{n-2}x_i}\mathrm dx_{n-1} \ldots \mathrm dx_1 \\= e^{-nx}\frac{x^{n-1}n^{n-2}}{(n-1)!}\,,
$$
where the volume integral may be evaluated by successive application of Leibniz's integral rule . 
A: One way to calculate $EM_n$ is to first calculate $P(M_n>m)$, for $m\in \mathbb R$, then use $EM_n=\int_0^\infty P(M_n>m)\,dm$.
Note that, in order for the event $\{M_n\le m\}$ to occur, it must be true that $X_1\le m$ and $X_2\le 2m-X_1$, and $X_3\le 3m-X_2-X_1$, and so on. Thus,
$$P(M_n\le m) = \int_0^m\int_0^{2m-x_1}\int_0^{3m-x_1-x_2}\cdots \int_0^{nm-x_1-\dots-x_{n-1}}e^{-(x_1+x_2+\cdots x_n)}dx_n\,dx_{n-1}\,\cdots dx_1\tag{*}$$
You can now use this expression to calculate $EM_n$. For example,
$$P(M_2\le m) = \int_0^me^{-x_1}\int_0^{2m-x_1}e^{-x_2}\,dx_2\,dx_1=\int_0^me^{-x_1}-e^{-2m}\,dx_1=1-me^{-2m}-e^{-m}$$
so, $$EM_2=\int_0^\infty P(M_2>m)\,dm=\int_0^\infty me^{-2m}+e^{-m}\,dm=\frac14+1$$
This does seem to suggest that $EM_n=\sum_1^n \frac1{k^2}$.
It might get messy, but perhaps you deduce a formula for $P(M_n\le m)$ by looking at more small cases, then prove it is correct by induction.
