# Covariance between an exponential random variable and the maximum of several exponential random variables

Suppose $X_1, \ldots, X_n$ are i.i.d. exponential RV with parameter $\lambda$. Let $Z = \max\{X_1, \ldots, X_6\}$. My goal is to find $\mathrm{cov}(X_1, Z)$.

I already know the c.d.f., p.d.f., and expectations of each $X_i$ and $Z$, I think I even found the joint p.d.f. of $Z$ and $X_1$ as:

$$f_{Z,\ X_1}(z,x_1) = \begin{cases} 5\lambda^2 e^{-\lambda x_1}(1-e^{-\lambda z})^4 & z \geq x_1 \geq 0 \\ 0 & \mathrm{otherwise} \end{cases}.$$ (This joint distribution is not correct, see below).

I found this by using $f_{Z,\ X_1} = f_{X_1} \ f_{Z|X_1}$ . But now I am sort of stuck because the integral $$E[X_1 Z] = \int\int_{R^2} x_1 z \ f_{Z,\ X_1} \ \mathrm{d}A$$ needed for $\mathrm{cov}(X_1,Z)$ gets cumbersome (and on my first attempt did not converge). Is there some slick thing to notice or do that I am missing? Note I am not yet confident with generating functions and moments.

• The total probability of your joint density $f_{Z,\ X_1}$ over the upper triangular region $z>x_1>0$ also diverges. Either it's just a typo or there's an error in your derivation that carried over to the calculation of $E[Z X_1]$. – Lee David Chung Lin May 31 '18 at 20:34
• The conditional distribution of $Z|X_1 = x_1$ is equivalent to $\max\{x_1, \max\{X_2, \ldots, X_6\}\}$. So it is a left-censored distribution, with a point mass on $x_1$. And this feature should appear in your joint as well. – BGM Jun 1 '18 at 3:01
• Thanks, here are my corrections to the joint distribution: $$$$\begin{split} F_{Z| \ X_1} &= P(Z \leq z \ | \ X_1 = x_1) \\ &= P(x_1 \leq X_1 \leq z \ | \ X_1 = x_1) \ldots P(x_1 \leq X_6 \leq z \ | \ X_1 = x_1) \\ &= \begin{cases} (1-e^{-\lambda(z-x_1)})^4 & z \geq x_1 \geq 0 \\ 0 & \mathrm{otherwise} \end{cases} \end{split}$$.$$ So with $f_{Z, \ X_1} = f_X \ f_{Z | \ X_1}$ I have: $$f_{Z, X_1} = \begin{cases} 5 \lambda^2 e^{-\lambda z} (1 - e^{-\lambda(z - x_1)})^4 & z \geq x_1 \geq 0 \\ 0 & \mathrm{otherwise} \end{cases}.$$ – frito_mosquito Jun 1 '18 at 14:27
• Even with the corrections my original question still stands, is there a better way to get to the eventual covariance? Are the better techniques confined to generating functions? – frito_mosquito Jun 1 '18 at 14:29
• @frito_mosquito (1) Your new derivation in the comment is still missing several things. In the short comment here I shall just say that the correct $f_{Z, X_1}$ involves the DiracDelta function and its 1st derivative. (2) Yes, there is a way to obtain the $E[Z X_1]$ without doing the integral with $f_{Z, X_1}$. The final result is $E[Z X_1] = \frac{343}{120 \lambda^2}$, where the coefficient $\approx 2.858333\ldots$. I will write up a post in a couple of days. – Lee David Chung Lin Jun 1 '18 at 19:58

The divergence of the integral as stated in the question is due to the presence of an atom (discrete point mass), as pointed out by BGM's comment. There are "proper" ways to deal with the point mass, but this answer will present the common technique that allows for an elementary solution.

Most of the analysis will be in the first section, and most of the calculation is left to the second. Two appendices come after the code block (numerical verification) to provide technical details.

## Strategy and Analysis

The analysis presented here applies to general $$n$$ and not just the requested $$n = 6$$.

Denote $$W$$ as the largest of the set of $$n-1$$ iid in contrast to $$Z$$ being the max of $$n$$. Here's let's say $$W$$ excludes $$X_1$$. $$W \equiv \max\{X_2,\,X_3,\,\ldots,\, X_6\}$$ Clearly $$W$$ is independent to $$X_1,\,$$ and its marginal density $$f_W$$ has the same form as $$f_Z$$ except with $$n-1$$ in place of $$n$$. The distribution of $$Z$$ as the max of $$n$$ iid exponential is well-known, as are the moments. Same for $$W$$.

Knowing that $$X_1$$ is interchangeable with any of the $$n-1$$ other $$X_i$$, we have the "$$\frac1n$$ and$$\frac{n-1}n$$" split that allows $$W$$ to come in. \begin{align*} E\bigl[ ZX_1] &= E\bigl[ZX_1~\big|\, Z = X_1 \bigr] \cdot\Pr\{Z = X_1\} + E\bigl[ZX_1~\big|\, Z > X_1\bigr] \cdot\Pr\{Z > X_1\} \\ &= E\left[Z^2~\middle|\, Z = X_1\right] \cdot \frac1n + E\bigl[ZX_1~\big|\, Z > X_1 \bigr] \cdot \frac{n-1}n \tag{1.a} \label{Eq_conditional_split} \\ &= E\left[Z^2 \right] \cdot \frac1n + E\bigl[WX_1~\big|\, W > X_1 \bigr] \cdot \frac{n-1}n \tag{1.b} \label{Eq_conditional_split_transfer} \end{align*} Going from the second line to the thrid line is trivial yet at the same time profound.

The first term $$E\left[Z^2~\middle|\, Z = X_1 \right] = E\left[Z^2 \right]$$ is effectively unconditional, because $$Z = X_i$$ for exactly one of the $$i = 1 \sim n$$ almost surely (ties has measure zero) while all $$X_i$$ play the same role. To be specific, due to $$E\left[Z^2~\middle|\, Z = X_i \right] = E\left[Z^2 \,\middle|\, Z = X_j \right]$$ for all $$i \neq j$$, we have \begin{align} E\left[Z^2\right] &= \sum_{i=1}^n \frac1n E\left[Z^2\,\middle|\, Z = X_i \right] \\ &= n \cdot \frac1n E\left[Z^2~\middle|\, Z = X_1 \right] = E\left[Z^2~\middle|\, Z = X_1 \right]\end{align} From the same line of reasoning, in fact we have identical distribution $$\left(Z^2~\middle|\, Z = X_1 \right) \overset{d}{=} Z^2$$. Note that this is conditioning on $$Z = X_1$$ and NOT conditioning on a given value of $$Z = X_1 = x_1$$.

As for the second term with $$X_1 < Z$$ going from Eq\eqref{Eq_conditional_split} to Eq\eqref{Eq_conditional_split_transfer}, by definition $$Z = \max\{W, X_1\}$$ so upon conditioning on $$Z > X_1$$ we have $$W = Z$$ as an identity between random variables (and not just identical in distribution).

Summary of the Strategy: diagonal split

The space can be split into disjoint cases as done in Eq\eqref{Eq_conditional_split}. For region $$Z = X_1$$ we bypass the joint density $$f_{X_1Z}$$ that is formidable by utilizing the marginal density $$f_Z$$ that is easy. For the other region, $$Z > X_1$$, we bypass the joint $$f_{X_1Z}$$ by recognizing its identical replacement $$f_{X_1W}$$, which is a joint density easily handled due to $$W \perp X_1$$.

This "diagonal split" applies to expectation of any well-behaving function $$g(X_1, Z)$$ in general and not just the product $$g(X_1, Z)=X_1 Z$$. Namely, we can compute $$E[g(X_1, Z)]$$ while not having to deal with $$f_{X_1Z}$$ which is cumbersome.

Nonetheless, $$\, f_{X_1Z}$$ is rather interesting and worth examining. Formally, we have the Heaviside step function (one form of indicator function), its derivative in the form of $$\delta(\cdot)$$ the Dirac delta function, and the first derivative of the Dirac delta function $$\delta'(\cdot)$$. Please consult wiki or Eq.(10) of the Wolfram page to see how $$\delta'(\cdot)$$ works. See Appendix.2 for an outline of the derivation of the full form of $$f_{X_1Z}$$ as in Eq\eqref{Eq_joint_density_full_form}.

## Calculation and Results

For the contribution from the "boundary" $$Z = X_1$$, one can obtain the moments with minor effort (see Appendix.1 for details if needed) \begin{align*} \text{At}~n = 6\,,& & E[Z] &= \frac{49}{20} \frac1{\lambda} & E\left[Z^2 \right] &= \frac{13489}{1800}\frac1{\lambda^2 } &&\qquad \tag{2} \label{Eq_Z_2nd_moment_n=6} \end{align*} For what it's worth, $$13489 = 7\cdot 41 \cdot 47$$ and $$13489/1800 \approx 7.49388888\ldots$$

In the $$Z > X_1$$ region, since $$W \perp X_1$$ the calculation of $$E\bigl[WX_1~\big|\, W > X_1 \bigr]$$ is easy to setup. The joint density is just $$f_{X_1W}(x, w) = f_X(x)\, f_W(w)$$ where the marginal density for $$X_1$$ is $$f_X$$ that is common to all $$X_i$$.

$$E\bigl[WX_1~\big|\, W > X_1 \bigr] = \frac{ E\bigl[W X_1\cdot \mathcal{I}_{W>X_1} \bigr] }{ \Pr\bigl\{ W > X_1 \bigr \} } = \frac{ \int_{x=0}^{\infty} \int_{w=x}^{\infty}~x \, w\, f_{X_1W}(x, w) \,\mathrm{d}w \,\mathrm{d}x }{ (n-1)/n }$$ where the $$\mathcal{I}_{\text{blah}}$$ is the indicator function, as in $$\mathcal{I}_{\text{blah}} = 1$$ wherever "blah" is true and $$\mathcal{I}_{\text{blah}} = 0$$ otherwise.

The denominator $$\Pr\{ W > X_1 \} = (n-1)/n$$ again due to symmetry (note that ties have zero probability marginally): $$\Pr\{ W > X_1 \} = 1 - \Pr\{ W \leq X_1 \} = 1 - \Pr\{ X_1~\text{is max} \} = 1 - \frac1n$$ For the record, below is the integral that one needs to evaluate. Note that the leading factor $$n/(n-1)$$ will cancel with that in Eq\eqref{Eq_conditional_split_transfer}, and this is not a coincidence. $$E\bigl[WX_1~\big|\, W > X_1 \bigr] = \frac{n}{n-1} \int_{x=0}^{\infty} \int_{w=x}^{\infty} \lambda^2 (n-1) x\,w\,e^{-\lambda(w+x)} \left(1 - e^{-\lambda w} \right)^{n-2} \,\mathrm{d}w \,\mathrm{d}x \tag{3.a} \label{Eq_E[XZ_on_X-smaller]_general_n}$$

The smaller $$n$$ cases are borderline manageable by hand, while in general one should not attempt this without CAS (computer algebra system) see the later half of Appendix.1. At $$n = 6$$, this evaluates to

$$E\bigl[ZX_1~\big|\, Z > X_1 \bigr] = \frac65 \frac{17381}{10800} \frac1{\lambda^2} \tag{3.b} \label{Eq_E[XZ_on_X-smaller]_n=6}$$ For what it's worth, $$17381 = 7 \cdot 13 \cdot 191$$ and $$17381/10800 \approx 1.609351851851851\ldots~$$, while $$6/5$$ times that is $$\approx 1.93122222\ldots$$

Take Eq\eqref{Eq_E[XZ_on_X-smaller]_n=6} and Eq\eqref{Eq_Z_2nd_moment_n=6} into Eq\eqref{Eq_conditional_split}, we have for $$n = 6$$:

\begin{align*} E[Z X_1] &= \frac16 E\left[ Z^2 \right] + \frac56 E\left[ Z X_1 ~\middle|\, X_1 < Z\right] \\ &= \frac{13489}{10800}\frac1{\lambda^2 } + \frac{17381}{10800} \frac1{\lambda^2}\\ &= \frac{343}{120} \frac1{\lambda^2} \\ &\approx \frac{ 2.8583333 }{\lambda^2} \end{align*} From here, one readily obtain the covariance using $$E[Z]$$ above and $$E[X_1] = 1 / \lambda$$. The covariance for $$n = 6$$ is $$49/(120 \lambda^2)$$. Using $$E\left[ Z^2 \right]$$ to get $$V[Z]$$ and do similar with $$E\left[ X_1^2 \right] = 1/\lambda^2$$, one obtains the correlation coefficient as $$\frac72\sqrt{7\,/\,767} \approx 0.3343639$$, moderately positive as expected.

Here's a quick and dirty Mathematica code verifying the results numerically (Monte Carlo).

ClearAll@SimulateEZX;
SimulateEZX[itrIn_: 17, InN_: 10^5, nIn_: 6, prcIn_: 20] := Module[{n (* # of iid *), itr, vN (* vector size *), itbyN, GTotal (* grand total*), data, Z, X, sumZX, prec, XisZ, sumBound, LenBound, bdN(* boundary case count *), bdTotal}, Timing[{itr, vN, n, prec} = Ceiling@{Norm@itrIn, Norm@InN, Norm@nIn, Norm@prcIn};  data = ConstantArray[0, {vN, n}]; Z = ConstantArray[0, vN];  X = Z; sumZX = ConstantArray[0, itr]; sumBound = sumZX;  LenBound = sumZX; XisZ = sumZX;
Do[ data = RandomVariate[ExponentialDistribution@1, {vN, n}, WorkingPrecision -> prec];
Z = Max /@ data; X = data[[;; , 1]];     sumZX[[k]] = Total[X Z];
XisZ[[k]] = Pick[X, MapThread[SameQ, {X, Z}]]; (* boundary case when X_1 \[Equal] Z *)
LenBound[[k]] = Length@XisZ[[k]];    sumBound[[k]] = Total[ XisZ[[k]]^2 ];
Print[{Row@{k, "th loop: "}, Row@{"sample_E[\!$$\*SubscriptBox[\(ZX$$, $$1$$]\)] = ", sumZX[[k]]/vN}, Row@{"  sample_E[\!$$\*SubscriptBox[\(ZX$$, $$1$$]\) | \ Z=\!$$\*SubscriptBox[\(X$$, $$1$$]\)] = ", sumBound[[k]]/LenBound[[k]]}}]
, {k, itr}];
{itbyN, GTotal, bdN, bdTotal} = {itr  vN, Total@sumZX, Total@LenBound, Total[sumBound]};
{ Row@{"Sample size = ", itbyN},  Row@{"  sample_E[\!$$\*SubscriptBox[\(ZX$$, $$1$$]\)] = ",      GTotal/itbyN} , Row@{"boundary case: count = ", bdN, " ,freq = ", N[bdN/itbyN, prec]}, Row@{"  sample_E[\!$$\*SubscriptBox[\(ZX$$, $$1$$]\) | \ Z=\!$$\*SubscriptBox[\(X$$, $$1$$]\)] = ", bdTotal/bdN},  Row@{"  sample_E[\!$$\*SubscriptBox[\(ZX$$, $$1$$]\) | \ \!$$\*SubscriptBox[\(X$$, $$1$$]\)<Z] = ", (GTotal - bdTotal)/(itbyN - bdN)}
} ]  ];
SimulateEZX[] (* default quick run  *)
SimulateEZX[20, 4.5 10^5, 6, 15] (* run with your own parameters *)


$$\color{white}.$$

## Appendix.1 : Distribution and Moments of Max of iid Exponential

Denote $$Z$$ as the maximum of $$n$$ random variables with iid exponential distribution with rate parameter $$\lambda$$. \begin{align*} F_Z(t) &= \left( 1 - e^{-\lambda t} \right)^n &&& f_Z(t) &= \lambda n e^{-\lambda t} \left( 1 - e^{-\lambda t} \right)^{n-1} \tag{A.1.1} \label{Eq_CDF_and_density_of_max} \end{align*} The moments can be obtained via routine (not necessarily easy) integration to be \begin{align*} E[Z] &= \frac1{\lambda} \mathcal{H}(n) & E\left[Z^2 \right] &= \frac1{ \lambda^2} \left( \frac{ \pi^2 }6 + \mathcal{H}(n)^2 -\psi'(n+1) \right) \tag{A.1.2} \label{Eq_Z_moments_general_n} \end{align*} where $$\mathcal{H}(n)\equiv \sum_k^n \frac1k$$ is the harmonic number and $$\psi'(\cdot)$$ is the first derivative of the digamma function (a.k.a polygamma with order 1).

A related calculation is the indispensable Eq\eqref{Eq_E[XZ_on_X-smaller]_general_n}. Ignoring the leading $$n/(n-1)$$ factor, this integral evaluates to a sequence $$S_n$$ that satisfies the following 2nd order difference equation

$$S_{n+1} = \frac{ n+ 1}{ (n-1) (n+2)^3 } \biggl[ n-4 + (n+2) \Bigl( 2 S_{n+1} \left(n^2+n-1 \right) - S_n + n^2 \Bigr) \biggr]$$

for $$n \geq 2$$, with the initial two terms $$S_2 = 1/2$$ (obvious) and $$S_3 = 47/54$$ (not so obvious). The first few terms are $$\{\frac12, \frac{ 47 }{ 54 }, \frac{ 335 }{ 288 }, \frac{ 12641 }{ 9000 }, \frac{ 17381 }{ 10800 }, \frac{ 1103219 }{ 617400 }, \ldots\}$$. Here for $$n = 6$$ what goes into Eq\eqref{Eq_E[XZ_on_X-smaller]_n=6} is $$S_6 = 17381/10800$$.

$$\color{white}.$$

## Appendix.2 : Derivation of Unconditional Joint Density $$f_{ZX_1}$$

I shall NOT include a sketch of the $$X_1$$-$$W$$ plane, with $$X_1$$ as the horizontal axis. Please make the sketch yourself as it will be helpful. The rectangular region that $$F_{X_1Z}(x,z)$$ covers is always tall-and-thin (height is as least as large as the width) but never short-and-fat.

Keep in mind that $$Z = \max\{W, X_1\}$$. From the sketch it should be easy to see that \begin{align*} F_{X_1Z}( x \, ,\, z ) &= \Pr\bigl\{ W \leq z ~~\& ~~X_1 \leq x \bigr\} &&\\ &= \Pr\bigl\{ W \leq z \bigr\} \cdot \Pr\bigl\{X_1 \leq x \bigr\} \qquad \because W \perp X_1 \\ &= F_W(z) \cdot F_X(x) \tag{A.2.1} \label{Eq_joint_CDF_elementary} \end{align*} where again $$F_W$$ is as seen in Eq\eqref{Eq_CDF_and_density_of_max} with $$W$$ instead of $$Z$$ and $$n-1$$ place of $$n$$.

Of course, $$F_{X_1Z}$$ should be different from $$F_{X_1W}$$, and one needs more explicit expressions.

$$F_{X_1Z}( x,\, z ) = F_W(z) \cdot F_X(x) \cdot \mathcal{I}_{z \geq x} \tag{A.2.2} \label{Eq_joint_CDF_indicator}$$

Even in its full indicator-form (with indicators regarding the range of $$X$$ and $$W$$), $$F_{X_1W} = F_W(z) \cdot F_X(x)$$ doesn't have this particular indicator $$\mathcal{I}_{z \geq x}$$.

Moving on. The indicator is a legitimate function, and the derivative of Eq\eqref{Eq_joint_CDF_indicator} follows the product rule as usual. \begin{align*} f_{X_1Z}( x \, ,\, z ) &= \frac{ \partial^2 F_W(z) F_X(x)}{ \partial x \partial z } \mathcal{I}_{z \geq x} + F_W(z) F_X(x)\frac{ \partial^2 \mathcal{I}_{z \geq x}}{ \partial x \partial z } \\ &\hspace{36pt} + \frac{ \partial F_W(z) F_X(x) }{ \partial z } \frac{ \partial \mathcal{I}_{z \geq x} }{ \partial x } + \frac{ \partial F_W(z) F_X(x) }{ \partial x } \frac{ \partial \mathcal{I}_{z \geq x} }{ \partial z } \end{align*} Apply the relations between the $$\{ \mathcal{I}_{\cdot},\delta(\cdot),\delta'(\cdot) \}$$ to eventually arrive at \begin{align*} f_{X_1Z}( x \, ,\, z ) &= \color{magenta}{ f_W(z)\, f_X(x)\, \mathcal{I}_{z \geq x} } + \delta(z-x) \biggl( \frac{ \partial }{ \partial z } \Bigl( F_W(z) F_X(x) \mathcal{I}_{z \geq x} \Bigr) \\ & \hspace{102pt} + F_W(z) f_X(x) - f_W(z) F_X(x) \biggr) \tag{A.2.3} \label{Eq_joint_density_full_form} \end{align*} The CDFs $$F_W, F_X$$ and their densities $$f_W$$, $$f_X$$ are established. One can pair up the densities in Eq\eqref{Eq_joint_density_full_form} to make it into a nice "matching" form if so desired.

The leading term (highlighted in magenta) corresponds to the region $$Z > X_1$$ that is identical to $$f_{X_1W}$$ conditionally.

• Very clever. If I may summarize a bit, the main ideas of your solution are to introduce the random variable $W$, use total expectation to split $E[ZX_1]$ into expectations conditioned on the value of $Z$, and then use independence of $W$ and $X_1$ to get a nice $f_{W,X_1}$ which can be integrated. – frito_mosquito Jun 4 '18 at 14:04
• @frito_mosquito Yes, that's the proper summary. And thanks for correcting my typo. – Lee David Chung Lin Jun 4 '18 at 14:15