pdf of $Z_{(n)}/Z_{(1)}$, where $Z \sim \text{Exponential}(1)$, i.i.d. (order statistics)

Question

Let $Z_i \sim \text{Exponential}(1)$, iid, $i=1,2, ..., n$ and let $Z_{(i)}$ be the $i^{\text{th}}$ order statistic.

How can I find the pdf of $T=Z_{(n)}/Z_{(1)}$?

(the ratio of maximum and minimum among $n$ random samples.)

---

The pdf of $(Z_{(1)},Z_{(n)})=(x,y)$ is

$f(x,y)=n(n-1)\ e^{-x}(e^{-x}-e^{-y})^{n-2} e^{-y} I(0<x<y)$,

and I could also find the pdf of $(Z_{(n)}/Z_{(1)},Z_{(1)})=(t,x)$, which is

$g(z,t)=n(n-1)\ x e^{-x}(e^{-x}-e^{-xt})^{n-2} e^{-xt} I(x>0,t>1)$.

What remains is to integrate g(z,t) by x, from 0 to inf. But I couldn't proceed.

I wish I could get some tips or answer.

Answer 1

Integrals of the form $\int xe^{-ax}$ can be solved using integration by parts. $$\int_{x_b}^{x_e} xe^{-ax} dx = \frac{-xe^{-ax}\Big|_{x_b}^{x_e} +\int_{x_b}^{x_e} e^{-ax}dx}{a} $$

and from $0$ to $\infty$ the first term on the right hand side cancels and the value is $\frac{1}{a^2}$

Your integral can be made in this form with $a$ the sum of your exponential terms $e^{c_1}e^{c_2}=e^{{c_1}+{c_2}}$. So you get as the integrand

$$x \sum \binom{n-2}{k} (-1)^k e^{-tx} \left(e^{-tx}\right)^{k}\left(e^{-x}\right)^{n-2-k}e^{-tx} $$

where we changed the (a+b)^n-2 product term in a sum

$$x \sum \binom{n-2}{k}(-1)^k e^{-(t+kt+(n-2-k)+t)x}$$

in which we changed the products like $e^{c_1}e^{c_2}$ into $e^{{c_1}+{c_2}}$ $$x \sum \binom{n-2}{k}(-1)^k e^{-(1+n-2+t+(t-1)k)x}$$

which rearranges the terms

and the integral of this integrand $\int_0^\infty$ is

$$\sum \binom{n-2}{k} \frac{(-1)^k}{(1+n-2+t+(t-1)k)^2}$$

which is of the form $\sum_{k=0}^{k=n} \binom{n}{k} \frac{(-1)^k}{(a+bk)^2}$ for which I do not easily find a similarity with other functions (the wolfram-alpha-site gives an expression with gamma functions) .

To check my calculus I did a verification by computation:

# settings
n <- 4
t <- 2
k <- c(0:n)

# numerical integration
integrate(function(x) x*exp(-x)*(exp(-x)-exp(-t*x))^n*exp(-t*x),lower=0,upper=Inf)

# integration by using the summation
sum(choose(n,k)*(-1)^k/(1+n+t+(t-1)*k)^2)

# using expression from wolfram alpha with gamma and digamma functions
a <- 1+n+t
b <- (t-1)
    # note: function is not universal and gamma of negative integers is not possible
(gamma(n+1)*gamma(a/b)*( digamma(a/b+n+1) - digamma(a/b) ) )/b^2/gamma(a/b+n+1)

which seems to work for all kinds of values $n$ and $t$

I've also checked your function, which is ok

# random data
sample1 <- matrix(rexp(10000*(n+2),1),10000)
# maximum and minimum
maxs <- apply(sample1,1,max)
mins <- apply(sample1,1,min)
# z_n/z_1
dif<-as.numeric(maxs/mins)
dif<-(dif<1000)*dif+(dif>=1000)*1000 #sensoring data >100

# plot historgram with calculated data
plot(hist(dif,breaks=seq(1,1000,1)),log="x")
lines(c(10:10000)/10,sapply(c(10:10000)/10,function(t) 10000*(n+2)*(n+1)*(t)^-0*sum(choose(n,k)*(-1)^k/(1+n+t+(t-1)*k)^2)),col=2)

Answer 2

A closed-form solution does appear to be possible, though perhaps the scenic route is needed.

Let $X \sim \text{Exponential}(1)$, let $X_1$ denote the sample minimum and $X_n$ denote the sample maximum in a sample of size $n$ drawn on parent $X$.

Then, it is straightforward (as the OP notes) that the joint pdf of $(X_1,X_n)$ is say $f(x_1, x_n)$:

where the Boole term here imposes the condition that $x_1 < x_n$.

Instead of seeking the pdf of of the desired ratio $\large\frac{X_n}{X_1}$ with domain of support on $(1,\infty)$, I am instead going to derive the pdf of the inverse $\large\frac{X_1}{X_n}$ which is bounded on (0,1) which seems to make the integration more manageable ... and then invert afterwards.

• Transforming $(X_1,X_n)$ to $(Y = \frac{X_1}{x_n}, Z = X_n)$ yields the joint pdf of $(Y,Z)$ as say $g(y,z)$:

where I am using the Transform function from the mathStatica package for Mathematica to automate, and domain[g] = {{y, 0, 1}, {z, 0, Infinity}}.

Then the marginal pdf of $Y$ is:

where:

• I add a cautionary note that there is a sign error in the solution returned (appears to be a Mathematica integration bug), which is corrected by replacing $(-1)^n$ with $(-1)^{n+1}$

• HarmonicNumber[n] denotes the $n^\text{th}$ harmonic number $H_n=\sum _{i=1}^n \frac{1}{i}$

Finally, making the transformation $R = \frac{1}{Y}$ yields the pdf of $R = \large \frac{X_n}{X_1}$ as:

---

$$\text{pdf}(r) = (-1)^{n+1} \frac{ n}{(r-1)^2 } \frac{ \Gamma (n) \Gamma \left(\frac{n r}{1-r}+1\right)}{\Gamma \left(\frac{n}{1-r}\right)} \left(H_{\frac{n r}{1-r}}-H_{\frac{n+r-1}{1-r}}\right) \quad \quad \text{for } r>1$$

---

Monte Carlo check

The following diagram compares:

• empirical pdf of Monte Carlo approximation of the max/min ratio when $n = 10$ (blue squiggly curve)
• the exact pdf derived above (when $n= 10$): dashed red curve underneath

All looks good.

Answer 3

For the joint pdf of order statistic we can use the following formula directly:

$g_{ij}(y_i,y_j)=\frac{n!}{(i-1)!(j-i-1)!(n-j)!}[F(y_i)]^{i-1}[F(y_j)-F(yi)]^{j-i-1}\times[1-F(y_j)]^{n-j}f(y_i)f(y_j) \tag {1}$

The poof can be found in Robert Hogg's book on Order statistics.$Y$ is the order statistics.

Your pdf for $Z_i$ is $f(z)=1*e^{-z}$. I suppose the support of your $pdf$ is $0<z_i<a$ which you should give for your question.

So $$F(x)=\int_{0}^xe^{-t}dt=1-e^{-x}$$ $$\therefore F(z_{(n)})=1-e^{-z_{(n)}}$$ $$ F(z_{(1)})=1-e^{-z_{(1)}}$$ $$ f(z_{(1)})=e^{-z_{(1)}}$$ $$ f(z_{(n)})=e^{-z_{(n)}}$$

Now let us plug in everything to $(1)$

$$f_{1,n}(z_{(1)},z_{(n)})=\frac{n!}{(n-2)!}[e^{-z_{(1)}}-e^{-z_{(n)}}]^{n-2}e^{-z_{(1)}}*e^{-z_{(n)}}\\=n(n-1)[e^{-z_{(1)}}-e^{-z_{(n)}}]^{n-2}e^{-z_{(1)}-z_{(n)}} \tag{2}$$

This is your joint $pdf$ for oder statistics $Z_{(1)}$ and $Z_{(n)}$

Next we do some variable tranform in oder to find pdf of $\frac{Z_{(n)}}{Z_{(1)}}$

Let $Y_1=\frac{Z_{(n)}}{Z_{(1)}}, Y_2=Z_{(1)}$ then $Z_{(n)}=Y_1Y_2, Z_{(1)}=Y_2$

The Jacobian is:

$$J=\begin{vmatrix} \frac{\partial z_{(n)}}{\partial y_1} & \frac{\partial z_{(n)}}{\partial y_2} \\ \frac{\partial z_{(1)}}{\partial y_1} & \frac{\partial z_{(1)}}{\partial y_2} \end{vmatrix}=y_2$$

Now the joint pdf of $Y_1$ and $Y_2$ can be found as: by plug in $(2)$

$$f_{Y_1,Y_2}(y_1,y_2)=f_{1n}(y_2,y_1 y_2)|J|=n(n-1)(e^{-y_2}-e^{-y_1y_2})^{n-2}e^{-y_2-y_1y_2}*y_2$$ Note we suppoes $0<z_i<a$ before so the $J=y_2$ is postive.

This is the joint pdf of $\frac{Z_{(n)}}{Z_{(1)}}=Y_1$ and $Z_{(1)}=Y_2$

Next,just integrate out $Y_2$ you will find the pdf of $\frac{Z_{(n)}}{Z_{(1)}}=Y_1$

$$f_{Y_1}(y_1)=\int_0^a f_{Y_1,Y_2}(y_1,y_2)dy_2\\=\int_0^an(n-1)(e^{-y_2}-e^{-y_1y_2})^{n-2}e^{-y_2-y_1y_2}*y_2dy_2\\=n(n-1)\int_0^a (e^{-y_2}-e^{-y_1y_2})^{n-2}*e^{-(1+y_1)y_2}y_2dy_2 \tag{3} $$

Ok, you are right, I cannot proceed from here, Martijn Weterings' method is a good one.