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

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.

• OK, the answer for the guys who will read this later: if you just want to calculate a robust scale estimator for a piece of data 1-install latest version of R 2-install the robustbase package 3-ready to go! but if you are developing a code outside this environment you need to use weighted high medians to minimize required calculations for Sn or Qn. – K-1 Oct 29 '10 at 7:28
• The link to the paper doesn't work. A proper reference (even better, with a quote of the most relevant information) would have helped us locate the information; as it stands it's useless when the link dies (as often happens). – Glen_b Mar 16 '16 at 1:49
• shouldn't it be k = h choose 2 = h(h-1)/2= 6 ? It doesn't change the end result, though. – a tiger Apr 25 '17 at 11:18
• why is Qn=Cn*2, why 2? how was it calculated? – lidox Feb 15 '18 at 15:12

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")

1. Load the package by typing:
library("robustbase")

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

• @user603: the tau you were referring to. Btw why isn't it widespread if it has such good statistical and computational efficiencies and breakdown point? – Quartz Nov 25 '13 at 13:01
• a) you can compute the mad and median online. From there it's trivial to compute the Tau. b) breakdown ain't robustness and the Tau has terrible bias in the presence of outliers. You can find more argument against it in section 5 of the Qn paper – user603 Nov 25 '13 at 13:26
• @user603 you mean this paper? wis.kuleuven.be/stat/robust/papers/publications-1994/… – German Demidov Mar 16 '16 at 0:05
• @user603 according to the paper, the bias curve tells us how much estimator can change due to a given fraction of contamination. $Q_n$ and $S_n$ were biased for my simulated examples (normal distribution + 20% of extremely high/low values), and the level of bias was comparable. May be I got something wrong, but both $S_n$ and $Q_n$ seem suffering from the same problem. – German Demidov Mar 16 '16 at 10:47
• @user603 sorry, the effect could not be seen for samples of size 100. I clearly see the problem using large sample sizes. All of them have terrible biases, but $\tau$ has the largest one. – German Demidov Mar 16 '16 at 10:59

(Very short answer) The text for commenting says

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, $$\{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).

• @Glen_b: go ahead and convert. Thank you for the clarification. – serv-inc Mar 16 '16 at 17:23
• @user603 The perhaps you could (as in the links in my comment) edit the essential information into the above answer -- as it stands at present it's not within the SE networks guidelines for answers. – Glen_b Mar 16 '16 at 23:17
• No problem, I will! (but it is really late here,) – user603 Mar 16 '16 at 23:51
• @user603 Thanks; I'll leave it here for now then – Glen_b Mar 17 '16 at 1:09

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={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;
}

• Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code. – gung - Reinstate Monica May 17 '18 at 2:02
• This is the naive O(n**2) algorithm~ – user603 Apr 6 '19 at 9:03