7.3 - joint distribution of two order statistics

Motivation

In the previous lecture we derived the marginal distribution of any \(j^{th}\) order statistic from a sample size \(n\) drawn i.i.d. from a population with pdf \(f_Y\) and CDF \(F_Y\).
We may often be interested in the joint distribution of two order statistics:
- To answer questions about joint probabilities
- To derive distribution of functions of two order statistics, e.g. \(R = Y_{(n)} - Y_{(1)}\)

Joint distribution of \(Y_{(i)}\), \(Y_{(j)}\)

For \(i<j\), we can heuristically derive \(f_{(i,j)}(x,w)\): “Probability of the \(i^{th}\) order statistic occurring at \(x\) and the \(j^{th}\) order statistic occurring at \(w\).”

Heuristic of the joint distribution of two order statistics

This leaves us with the following multinomial “buckets”:

multinomial buckets for the joint distribution of the i and jth order statistics

Thus the joint pdf of \((Y_{(i)},Y_{(j)})\) for \(i<j\) follows:

\[ f_{(i,j)}(x,w) = \frac{n!}{(i-1)!1!(j-i-1)!1!(n-j)!}F_Y(x)^{i-1}f_Y(x) \left(F_Y(w)-F_Y(x)\right)^{j-i-1} f_Y(w) \left(1-F_Y(w)\right)^{n-j}\]

\[ = \frac{n!}{(i-1)!(j-i-1)!(n-j)!}F_Y(x)^{i-1}f_Y(x) \left(F_Y(w)-F_Y(x)\right)^{j-i-1} f_Y(w) \left(1-F_Y(w)\right)^{n-j},x,w\in population\ support\ and\ x<w\]

Example: Joint distribution of exponential min and max

Let \(Y_1\), \(Y_2\),…,\(Y_n\) be an i.i.d. sample from an \(EXP(\lambda)\) population.
Thus:

\[F_Y(y) = 1-e^{-\lambda y}\]

\[f_Y(y) = \lambda e^{-\lambda y }\]

Find joint pdf of \(Y_{(1)}\) and \(Y_{(n)}\).

\[\scriptsize f_{(1,n)}(x,w) = \frac{n!}{(1-1)!(n-1-1)!(n-n)!}F_Y(x)^{1-1}f_Y(x) \left(F_Y(w)-F_Y(x)\right)^{n-1-1} f_Y(w) \left(1-F_Y(w)\right)^{n-n}\]

\[\scriptsize =n(n-1)\lambda e^{-\lambda x} \left(e^{- \lambda x} - e^{-\lambda w}\right)^{n-2}\lambda e^{-\lambda w} = n(n-1)\lambda^2 e^{-\lambda (x+w)} \left(e^{- \lambda x} - e^{-\lambda w}\right)^{n-2}, 0<x<w<\infty\]

Application: find pdf of \(\small R = Y_{(n)} - Y_{(1)}\).

Have: \(f_{(1,n)}(x,w)\), joint pdf of \((Y_{(1)}, Y_{(n)})\)
Need: pdf of a transformation \(R = Y_{(n)} - Y_{(1)}\)
Approach: find joint pdf of \(2\rightarrow 2\) transformation \((Y_{(1)}, Y_{(n)}) \rightarrow (R,S)\) for some \(S\), integrate out \(S\) to find marginal of \(R\).
Required method: Jacobian method for \(2 \rightarrow 2\) transformations

Refresher: Jacobian method for \(2 \rightarrow 2\) transformations

Let \((X,Y)\) be jointly continuous with joint pdf \(f_{X,Y}(x,y)\)
Define \(2 \rightarrow 2\) transformation: \(U = g(X,Y), V = h(X,Y)\)
Then:

\[f_{U,V}(u,v) = f_{X,Y}(g^{-1}(u,v), h^{-1}(u,v)) |det(J)|,\]

where:

\[ \scriptsize J = \begin{bmatrix} \frac{\partial g^{-1}}{\partial u} & \frac{\partial g^{-1}}{\partial v} \\ \frac{\partial h^{-1}}{\partial u} & \frac{\partial h^{-1}}{\partial v} \end{bmatrix}\]

\[\scriptsize |det(J)| = \left|\frac{\partial g^{-1}}{\partial u}\frac{\partial h^{-1}}{\partial v} - \frac{\partial h^{-1}}{\partial u} \frac{\partial g^{-1}}{\partial v}\right|\]

Specifying the \(2 \rightarrow 2\) transformation

Let \(R = g(Y_{(1)},Y_{(n)}) = Y_{(n)} - Y_{(1)}\), \(S=h(Y_{(1)},Y_{(n)})=Y_{(1)}\)
Find joint pdf \(f_{R,S}(r,s)\).

Joint support of \((Y_{(1)},Y_{(n)})\): \(0 < Y_{(1)}<Y_{(n)} < \infty\)
Joint support mapping:

\(Y_{(1)} = 0 \Rightarrow R = Y_{(n)}, S = 0\)
\(Y_{(n)} = Y_{(1)} \Rightarrow R = 0, S = Y_{(1)}\)

\[\therefore R > 0, S> 0\]

Inverses and Jacobian

\[R = g(Y_{(1)},Y_{(n)}) = Y_{(n)} - Y_{(1)}\] \[S=h(Y_{(1)},Y_{(n)})=Y_{(1)}\]

\[\Rightarrow Y_{(1)} = g^{-1}(R,S) = S\] \[\Rightarrow Y_{(n)} = h^{-1}(R,S) = R+S\]

\[ J = \begin{bmatrix} \frac{\partial g^{-1}}{\partial r} & \frac{\partial g^{-1}}{\partial s} \\ \frac{\partial h^{-1}}{\partial r} & \frac{\partial h^{-1}}{\partial s} \end{bmatrix}= \begin{bmatrix} 0 & 1 \\ 1 &1 \end{bmatrix}\]

\[ |det(J)| = |-1| = 1\]

Plugging and chugging:

\[f_{R,S}(r,s) = f_{(1,n)}(s,r+s)\cdot 1\]

\[ = n(n-1)\lambda^2 e^{-\lambda (s+(r+s))} \left(e^{- \lambda s} - e^{-\lambda (r+s)}\right)^{n-2}\]

\[ = n(n-1)\lambda^2 e^{-2\lambda s} e^{-\lambda r} \left(e^{- \lambda s}\left(1- e^{-\lambda r}\right)\right)^{n-2}, r>0, s>0\]

Finding marginal of \(R\)

Recall our goal is just to find \(f_R(r)\):

\[ f_R(r)= \int_0^\infty n(n-1)\lambda^2 e^{-2\lambda s} e^{-\lambda r} \left(e^{- \lambda s}\left(1- e^{-\lambda r}\right)\right)^{n-2}\, ds\]

\[ =n(n-1)\lambda^2e^{-\lambda r} \left(1- e^{-\lambda r}\right)^{n-2} \int_0^\infty e^{-2\lambda s} \left(e^{- \lambda s}\right)^{n-2}\, ds\]

\[ =n(n-1)\lambda^2e^{-\lambda r} \left(1- e^{-\lambda r}\right)^{n-2} \int_0^\infty e^{-2\lambda s -n\lambda s+2\lambda s} \, ds\]

\[=n(n-1)\lambda^2e^{-\lambda r} \left(1- e^{-\lambda r}\right)^{n-2} \underbrace{\int_0^\infty e^{ -n\lambda s} \, ds}_{kernel\ of\ EXP(n \lambda)\ integrated}\]

\[=n(n-1)\lambda^2e^{-\lambda r} \left(1- e^{-\lambda r}\right)^{n-2}\cdot \frac{1}{n\lambda} = (n-1)\lambda e^{-\lambda r} \left(1- e^{-\lambda r}\right)^{n-2}, r > 0\]

Verifying via simulation study

library(purrrfect)
library(tidyverse)



(exponential_range_sim <- parameters(~n, ~lambda,
                             c(2, 5, 10), c(0.5, 1.0, 1.5)
                            )
                  %>%  add_trials(10000)
                  %>%  mutate(y_sample = pmap(list(n, lambda), .f = \(n, l) rexp(n, rate = l)))
                  %>%  mutate(y1 = map_dbl(y_sample, .f = min),
                              yn = map_dbl(y_sample, .f = max),
                              r = yn - y1)
                  %>%  mutate(f_r =  (n-1)*lambda*exp(-lambda *r) *(1-exp(-lambda*r))^(n-2) 
                              )
)  %>% head()

# A tibble: 6 × 8
      n lambda .trial y_sample      y1    yn     r    f_r
  <dbl>  <dbl>  <dbl> <list>     <dbl> <dbl> <dbl>  <dbl>
1     2    0.5      1 <dbl [2]> 0.577   3.87 3.29  0.0964
2     2    0.5      2 <dbl [2]> 0.138   4.45 4.31  0.0580
3     2    0.5      3 <dbl [2]> 0.623   2.60 1.97  0.186 
4     2    0.5      4 <dbl [2]> 1.10    3.29 2.19  0.167 
5     2    0.5      5 <dbl [2]> 0.0269  2.45 2.42  0.149 
6     2    0.5      6 <dbl [2]> 1.32    1.76 0.441 0.401

Plotting results

ggplot(data = exponential_range_sim, aes(x = r))+
  geom_histogram(aes(y = after_stat(density)), 
                 fill = 'goldenrod',
                 center = 0.1, binwidth = .2) +
  geom_line(aes(y = f_r), col='cornflowerblue') +
  facet_grid(n~lambda, labeller = label_both, scales='free_y')+
  theme_classic(base_size = 16) +
  labs(title='Simulated and analytic densities for range of exponentials')

simulated and analytic densities of range of exponentials

Example: joint distribution of uniform order stats

Let \(Y_1\), \(Y_2\),…,\(Y_5\) be an i.i.d. sample of \(UNIF(0,1)\) random variables
Find \(P(Y_{(2)} < 0.3, Y_{(4)} > 0.7)\).
Parent population CDF and pdf, for \(0 < y < 1\):

\[F_{Y}(y)= y\] \[f_{Y}(y)= 1\]

\[f_{(2,4)}(x,w) = \frac{5!}{(2-1)!(4-2-1)!(5-4)!}x^1\cdot 1 \cdot (w-x)^1\cdot 1\cdot (1-w)^1\]

\[=120x(w-x)(1-w), 0 < x < w < 1\]

Finding \(P(Y_{(2)} < 0.3, Y_{(4)} > 0.7)\)

Region of integration:

\[P(Y_{(2)} < 0.3, Y_{(4)} > 0.7) = \int_{0.7}^1\int_0^{0.3} 120x(w-x)(1-w)\, dx \, dw\]

\[ = \int_{0.7}^1(1-w)(60x^2w - 40x^3)\big|_{x=0}^{0.3}\, dw\]

\[ = \int_{0.7}^1(1-w)(5.4w-1.08)\, dw = 0.1458\]

Simulating

Notes:
- As we do not have a grid to simulate over (\(n = 5\), and parent population remains \(UNIF(0,1)\) for each trial), we can simply use replicate() from the purrrfect package.
- For general order statistics, we must use map() to create a list column of the sorted sample.
- map_dbl(list_column, pluck(j)) and map_dbl(list_column, j) are identical, simply “plucking” the \(j^{th}\) element from each vector in the list.

(uniform_order_stats <- replicate(10000,  
                                  runif(5, min = 0, max = 1),
                                  .as = y_sample
                                  )
                  %>% mutate(y_sorted = map(y_sample, sort))
                  %>% mutate(y2 = map_dbl(y_sorted, pluck(2)),
                             y4 = map_dbl(y_sorted, 4)
                             )
)  %>% head()

# A tibble: 6 × 5
  .trial y_sample  y_sorted     y2    y4
   <dbl> <list>    <list>    <dbl> <dbl>
1      1 <dbl [5]> <dbl [5]> 0.413 0.601
2      2 <dbl [5]> <dbl [5]> 0.341 0.711
3      3 <dbl [5]> <dbl [5]> 0.341 0.620
4      4 <dbl [5]> <dbl [5]> 0.620 0.831
5      5 <dbl [5]> <dbl [5]> 0.309 0.361
6      6 <dbl [5]> <dbl [5]> 0.545 0.826

Aggregating to find simulated probability

(uniform_order_stats
  %>% summarize('P(Y2 < 0.3, Y4 > 0.7)' = mean(y2 < 0.3 & y4 > 0.7))
)

# A tibble: 1 × 1
  `P(Y2 < 0.3, Y4 > 0.7)`
                    <dbl>
1                   0.145