Jak obliczyć estymator skali Qn Rousseeuw'a i Crouxa (1993) dla dużych próbek?

13

Niech więc dla bardzo krótkiej próbki, takiej jak , można ją obliczyć od znalezienia statycznego tego rzędu różnicy par: { 1 , 3 , 6 , 2 , 7 , 5 } kQn=Cn.{|XiXj|;i<j}(k){1,3,6,2,7,5}k

    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

ZatemQn=Cn.2

Oczywiście w przypadku dużych próbek, które mówią, że zawierają 80 000 rekordów, potrzebujemy bardzo dużej pamięci.

Czy w ogóle można obliczyć w przestrzeni 1D zamiast 2D?Qn

Link do odpowiedzi ftp://ftp.win.ua.ac.be/pub/preprints/92/Timeff92.pdf, chociaż nie rozumiem tego w pełni.

K-1
źródło
1
OK, odpowiedź dla facetów, którzy przeczytają to później: jeśli chcesz tylko obliczyć solidny estymator skali dla kawałka danych 1-zainstaluj najnowszą wersję R 2-zainstaluj pakiet robustbase 3-gotowy! ale jeśli tworzysz kod poza tym środowiskiem, musisz użyć ważonych wysokich median, aby zminimalizować wymagane obliczenia dla Sn lub Qn.
K-1
1
Link do artykułu nie działa. Właściwe odniesienie (nawet lepsze, z podaniem najbardziej istotnych informacji) pomogłoby nam je znaleźć; w obecnej postaci jest bezużyteczny, gdy umiera link (jak to często bywa).
Glen_b
2
czy nie powinno być k = h wybrać 2 = h (h-1) / 2 = 6 ? Nie zmienia to jednak wyniku końcowego.
tygrys
dlaczego Qn = Cn * 2, dlaczego 2? jak to zostało obliczone?
lidox

Odpowiedzi:

15

Aktualizacja: sedno problemu polega na tym, że aby osiągnąć złożoność czasową O(nlog(n)) , potrzeba kolejności przechowywania O(n) .


Nie, O(nlog(n)) jest dolną teoretyczną granicą złożoności czasowej (patrz (1)) wyboru elementu kth spośród wszystkich n(n1)2 możliwe|xixj|:1i<jn.

Można dostać O(1) miejsca, ale tylko przez naiwnie sprawdzając wszystkie kombinacje xixj w czasie O(n2) .

Dobrą wiadomością jest to, że można użyć estymatora skali τ (patrz (2) i (3) dla ulepszonej wersji i niektórych porównań czasowych), zaimplementowanego w funkcji scaleTau2()w Rpakiecie robustbase. Jednoznaczny estymator τ jest dwuetapowym (tj. Ponownie ważonym) estymatorem skali. Ma 95 procent wydajności Gaussa, 50 procent punktu przebicia oraz złożoność czasu O(n) i przestrzeni O(1) (a ponadto można go łatwo ustawić w trybie online, zmniejszając o połowę koszty obliczeniowe przy wielokrotnym użyciu - chociaż ty będzie musiał zagłębić się w Rkod, aby wdrożyć tę opcję, jest to raczej proste do zrobienia).

  1. Złożoność wyboru i rankingu w X + Y i macierzach z posortowanymi kolumnami GN Frederickson i DB Johnson, Journal of Computer and System Sciences Tom 24, Issue 2, April 1982, Pages 197-208.
  2. Yohai, V. i Zamar, R. (1988). Oszacowania regresji dla wysokiego punktu awarii za pomocą minimalizacji skutecznej skali. Journal of American Statistics Association 83 406–413.
  3. Maronna, R. i Zamar, R. (2002). Solidne oszacowania lokalizacji i dyspersji dla zbiorów danych wielowymiarowych. Technometria 44 307–317

Edytuj Aby użyć tego

  1. Odpal się R(jest bezpłatny i można go pobrać tutaj )
  2. Zainstaluj pakiet, wpisując:
install.packages("robustbase")
  1. Załaduj pakiet, wpisując:
library("robustbase")
  1. Załaduj plik danych i uruchom funkcję:
mydatavector <- read.table("address to my file in text format", header=T)
scaleTau2(mydatavector)
użytkownik603
źródło
2
@ user603: tau, o którym mówiłeś. Przy okazji, dlaczego nie jest szeroko rozpowszechniony, jeśli ma tak dobrą wydajność statystyczną i obliczeniową oraz punkt podziału?
Kwarc
2
a) możesz obliczyć szaloną i medianę online . Stamtąd obliczenie Tau jest banalne. b) rozpad nie jest odporny, a Tau ma straszne nastawienie w obecności wartości odstających. Możesz znaleźć więcej argumentów przeciwko temu w sekcji 5 artykułu Qn
603
1
@ user603 masz na myśli ten papier? wis.kuleuven.be/stat/robust/papers/publications-1994/…
Niemiecki Demidov
1
QnSnSnQn
1
τ
0

(Bardzo krótka odpowiedź) Tekst do komentowania mówi

unikaj odpowiadania na pytania w komentarzach.

Qn

EDYTOWAĆ

Qn

{xi}i=1Nn<N{xi}i=tn+1tQnNn+1Qn{Qni}i=1Nn+1

Qni|Qni1 O(nlog(n))Qni

Qn{xi}i=1NO(n2)

serv-inc
źródło
Chociaż nie powinieneś odpowiadać w komentarzach, nie powinieneś również zamieszczać komentarzy jako odpowiedzi, a jeśli twoja odpowiedź jest tylko linkiem, nie jest to odpowiedź (ale może być komentarzem). Jeśli chcesz, aby była to odpowiedź, a nie komentarz, twoja odpowiedź powinna w jakiś sposób zawierać istotne informacje, takie jak cytat z prawidłowo przywołanego linku lub własne wyjaśnienie ważnych szczegółów. Jeśli możesz, podaj niezbędne dane; alternatywnie mogę przekonwertować to na komentarz dla ciebie.
Glen_b
@Glen_b: śmiało i konwertuj. Dziękuję za wyjaśnienie.
serv-inc
1
@ user603 Być może mógłbyś (tak jak w linkach w moim komentarzu) edytować istotne informacje w powyższej odpowiedzi - w obecnym stanie nie jest to zgodne z wytycznymi sieci SE dotyczącymi odpowiedzi.
Glen_b
Nie ma problemu, zrobię to! (ale tutaj jest naprawdę późno)
użytkownik603
@ user603 Dzięki; Na razie zostawię to tutaj
Glen_b
0

to moja implementacja Qn ...

Programowałem to w C, a wynik jest następujący:

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;
}
zwycięzca
źródło
1
Chociaż implementacja jest często pomieszana z treścią merytoryczną w pytaniach, mamy być witryną do dostarczania informacji o statystykach, uczeniu maszynowym itp., A nie o kodzie. Dobrze jest również podać kod, ale proszę opracować merytoryczną odpowiedź w tekście dla osób, które nie czytają tego języka wystarczająco dobrze, aby rozpoznać i wyodrębnić odpowiedź z kodu.
Gung - Przywróć Monikę
Jest to naiwny algorytm O (n ** 2) ~
użytkownik603