How to calculate Rousseeuw’s and Croux’ (1993) Qn scale estimator for large samples?

Question

Let $Q_n = C_n.\{|X_i-X_j|;i < j\}_{(k)}$ so for a very short sample like $\{1,3,6,2,7,5\}$ it can be calculated from finding the $k$th order static of pairwise differences:

    7 6 5 3 2 1
1   6 5 4 2 1
2   5 4 3 1
3   4 3 2
5   2 1
6   1
7

h=[n/2]+1=4

k=h(h-1)/2=8

Thus $Q_n=C_n. 2$

Obviously for large samples saying consist of 80,000 records we need very large memory.

Is there anyway to calculate $Q_n$ in 1D space instead of 2D?

A link to the answer ftp://ftp.win.ua.ac.be/pub/preprints/92/Timeff92.pdf although I cannot fully understand it.

Answer 1

Update: The crux of the problem is that in order to achieve the $O(n\log(n))$ time complexity, one needs in the order of $O(n)$ storage.

---

No, $O(n\log(n))$ is the lower theoretical bound for the time complexity of (see (1)) selecting the $k^{th}$ element among all $\frac{n(n-1)}{2}$ possible $|x_i - x_j|: 1 \leq i \lt j \leq n$.

You can get $O(1)$ space, but only by naively checking all combinations of $x_i-x_j$ in time $O(n^2)$.

The good news is that you can use the $\tau$ estimator of scale (see (2) and (3) for an improved version and some timing comparisons), implemented in the function scaleTau2() in the R package robustbase. The univariate $\tau$ estimator is a two-step (i.e. re-weighted) estimator of scale. It has 95 percent Gaussian efficiency, 50 percent breakdown point, and complexity of $O(n)$ time and $O(1)$ space (plus it can easily be made 'online', shaving off half the computational costs in repeated use -- although you will have to dig into the R code to implement this option, it is rather straightforward to do).

1. The complexity of selection and ranking in X + Y and matrices with sorted columns G. N. Frederickson and D. B. Johnson, Journal of Computer and System Sciences Volume 24, Issue 2, April 1982, Pages 197-208.
2. Yohai, V. and Zamar, R. (1988). High breakdown point estimates of regression by means of the minimization of an efficient scale. Journal of the American Statistical Association 83 406–413.
3. Maronna, R. and Zamar, R. (2002). Robust estimates of location and dispersion for high- dimensional data sets. Technometrics 44 307–317

Edit To use this

1. Fire up R (it's free and can be downloaded from here)
2. Install the package by typing:

install.packages("robustbase")

3. Load the package by typing:

library("robustbase")

4. Load your data file and run the function:

mydatavector <- read.table("address to my file in text format", header=T)
scaleTau2(mydatavector)

Answer 2

(Very short answer) The text for commenting says

avoid answering questions in comments.

so here it goes: There is a paper about an online algorithm which seemingly runs quite well: Applying the $Q_n$ Estimator Online.

EDIT

(by user user603). The algorithm linked to in this article is a moving window version version of the $Q_n$.

Given a large sample $\{x_i\}_{i=1}^N$ divided into time windows of width $n<N$, $\{x_i\}_{i=t-n+1}^t$ we can apply the $Q_n$ to each time window yielding $N-n+1$ values of the $Q_n$. Denote these values $\{Q_n^i\}_{i=1}^{N-n+1}$

The algorithm cited here allows to obtain $Q_n^i|Q_n^{i-1}$ at an average cost less than the worst case $O(n\log(n))$ needed to compute $Q_n^i$ from scratch.

This algorithm can however not be used to compute the $Q_n$ of the full original sample $\{x_i\}_{i=1}^N$. It also needs to maintain an buffer whose size can be as large as $O(n^2)$ (though it is often much smaller).

Answer 3

this is my implement of Qn...

I was programming this in C and the result is this:

void bubbleSort(double *datos, int N)
{
 for (int j=0; j<N-1 ;j++)     
  for (int i=j+1; i<N; i++)    
   if (datos[i]<datos[j])      
   {
    double tmp=datos[i];
    datos[i]=datos[j];
    datos[j]=tmp;
   }
}

double  fFactorial(long N)    
{
 double factorial=1.0;

 for (long i=1; i<=N; ++i)
  factorial*=(double)i;

 return factorial;  
}

double fQ_n(double *datos, int N)  // Rousseeuw's and Croux (1993) Qn scale estimator
{
 bubbleSort(datos, N);

 int m=(int)((fFactorial((long)N))/(fFactorial(2)*fFactorial((long)N-2)));

 double D[m];
 //double Cn=2.2219;      //not used now :) constant value https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/qn_scale.htm

 int k=(int)((fFactorial((long)N/2+1))/(fFactorial(2)*fFactorial((long)N/2+1-2)));

 int y=0;

 for (int i=0; i<N; i++)
  for (int j=N-1; j>=0; j--)
   if (i<j)
   {
    D[y]=abs(datos[i]-datos[j]);
    y++;
   }

 bubbleSort(D, m);

 return D[k-1];
}

int main(int argc, char **argv)    
{
 double datos[6]={1,2,3,5,6,7};
 int N=6;

 // Priting in terminal the final solution
 printf("\n==[Results] ========================================\n\n");

 printf(" Q_n=%0.3f\n",fQ_n(datos,N));

 return 0;
}