Order-Stat Talk

I remember of one of my third year courses, Order Statistics and Non-Parametric methods. One day we were told about an assignment which was really fun attempting! I was just thinking if I can write and share about my experience of doing that task with my very little knowledge of statistics and programming. Happy reading :)!

As the statisticians are always in need of the tables of many particular distributions for further mathematical manipulation, which may not be very available, they can at least attempt to produce them by themselves. But, I faced some of silly troubles while doing this work. I hope to get advice about them from the respected statisticians roaming around :). Nonetheless, I was happy with my job at the end, but couldn’t find an easier way to make it done.

Before I move on in my article, let me first tell you about the task I am talking about:
Consider $X_1, X_2, \ldots, X_n$ be a random sample taken from standard exponential distribution. Here, $n = 1,\ldots,5.$ Find the means and variances of ith order statistics $X_{i:n}.$

A brief introduction to Order Statistics:

Let, $X_1, X_2,\ldots,X_n$ be a random sample from a probability density function $f(x)$ and corresponding cumulative distribution function be denoted by $F(x).$ Let us consider the sorting of the sample observations such that,

$X_{i:n} \le X_{2:n} \le \ldots \le X_{n:n}$

Here $X_{i:n}$ is known as the ith Order Statistic in a sample of size $n.$

Distribution of the single Order Statistics:

$$f_{i:n} = \frac{n!}{(i-1)!(n-i)!} [F(x)]^{i-1} [1-F(x)]^{n-i} f(x), \quad -\infty < x< \infty. $$

Following this we get $f_{i:n} (x)$ for standard exponential distribution as:

$$f_{i:n}(x) = \frac{n!}{(i-1)!(n-i)!}[1- e^{x_{i:n}}]^{i-1}[e^{-x_{i:n}(n-i+1)}], \quad x>0$$

Here comes the programming part. First save the following commands in an external file, raw1.txt, say in drive e:

fnc<-function(p){for(n in 1:5){
for(i in 1:n){
a<-(factorial(n)/(factorial(i-1)*factorial(n-i)))
integrand<-function(x) {a*(x^p)*((1-exp(-x))^(i-1))*exp(-x*(n-i+1))}
r<-integrate(integrand, lower=0, upper=Inf)
print(r)}}
f1<-fnc(1)
print(f1)
f2<-fnc(2)
print(f2)

Using the following commands we can save the output (the first raw moments or the means) in a text file out1.txt:

sink(“e:\\out1.txt”)
source(“e:\\raw1.txt”)
sink( )

The out1.txt file contains:

1 with absolute error < 6.4e-06
0.5 with absolute error < 8.6e-06
1.5 with absolute error < 1.3e-05
0.3333333 with absolute error < 8e-08
0.8333333 with absolute error < 2.7e-05
1.833333 with absolute error < 1.9e-05
0.25 with absolute error < 3.2e-05
0.5833333 with absolute error < 3.5e-07
1.083333 with absolute error < 5.5e-05
2.083333 with absolute error < 2.5e-05
0.2 with absolute error < 8.3e-05
0.45 with absolute error < 2.9e-08
0.7833333 with absolute error < 9.5e-07
1.283333 with absolute error < 9.6e-05
2.283333 with absolute error < 3.2e-05
NULL
2 with absolute error < 7.1e-05
0.5 with absolute error < 4e-05
3.5 with absolute error < 0.00014
0.2222222 with absolute error < 1.5e-06
1.055556 with absolute error < 0.00013
4.722222 with absolute error < 0.00021
0.125 with absolute error < 4.1e-05
0.5138889 with absolute error < 6.4e-06
1.597222 with absolute error < 2.4e-07
5.763889 with absolute error < 0.00028
0.08 with absolute error < 8.2e-06
0.305 with absolute error < 4.1e-07
0.8272222 with absolute error < 1.7e-05
2.110556 with absolute error < 4e-07
6.677222 with absolute error < 0.00035


The values are actually the first and second raw moments of the order statistics, given with the absolute errors. This list of means can not be read as a numeric vector in R. So we have to modify them by deleting the absolute errors manually. Also, we need to make two distinct external files for the first and second raw moments for using them further. I found it really problematic. This can not be a solution for a large value of n! Need suggestions on this issue!

Now after reading those moments back in R, we need to calculate the variance by the rule of moments:

$$Var(x)= E(x^2) - {E(x)}^2 $$

The lastly produced table will be something like this:

Table: Means and Variances of Exponential Order Statistics for $n$ up to 5

n i mu_{i:n}sigma_{i, i:n}
1 1 1 1.0000000
2 1 0.5 0.2500000
2 2 1.5 1.2500000
3 1 0.333 0.1111111
3 2 0.833 0.3611116
3 3 1.833 1.3611121
4 1 0.25 0.0625000
4 2 0.583 0.1736112
4 3 1.083 0.4236116
4 4 2.083 1.4236126
5 1 0.2 0.0400000
5 2 0.45 0.1025000
5 3 0.783 0.2136111
5 4 1.283 0.4636124
5 5 2.283 1.4636124


Reference: A First Corse in Order Statistics, Barry C Arnold, N. Balakrishnan, H.N. Nagaraja, John Wiley & Sons.