4
$\begingroup$

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.

$\endgroup$
1
4
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ It would be neater if you could post the result/plot of your simulation (compared to theoretical derivation), rather than the code to generate it $\endgroup$ – wolfies Oct 5 '17 at 5:13
  • $\begingroup$ The code is more like an appendix for a sanity check and not really a result. I placed it such that anyone can do the check as well. It is a rather simple copy and past. I did this code to test my result as well as the second block of code to test the g(z,t) equation by 456123. ----- The output is boring. The first block is just an output of three times the same numbers. The second block is a histogram (from a random sample) with a density function on top of it (modeled). Neither very special. $\endgroup$ – Sextus Empiricus Oct 5 '17 at 11:23
5
$\begingroup$

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)$:

enter image description here

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)$:

enter image description here

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:

enter image description here

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

enter image description here All looks good.

$\endgroup$
1
  • $\begingroup$ The harmonic number is in this case a generalized harmonic number, for real values. $\endgroup$ – Sextus Empiricus Oct 5 '17 at 18:12
0
$\begingroup$

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.

$\endgroup$
10
  • 1
    $\begingroup$ What is parameter "a"? Why should an Exponential distribution have an upper finite bound? $\endgroup$ – wolfies Oct 5 '17 at 4:07
  • 1
    $\begingroup$ I see, I made a mistake, he means $e^{-x}$ I will edit my answer $\endgroup$ – Deep North Oct 5 '17 at 4:19
  • 1
    $\begingroup$ I added that part. I suppose $0<z_i<a$. OP did not tell us the domain of $pdf$. I made that part. OP should tell us in which domain this pdf hold. $\endgroup$ – Deep North Oct 5 '17 at 5:37
  • 1
    $\begingroup$ For oder statistic there should have a upper limit there. $\endgroup$ – Deep North Oct 5 '17 at 5:45
  • $\begingroup$ I think the third line and fourth line of equation (3) doesn't match. $\endgroup$ – 456 123 Oct 5 '17 at 7:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.