Medcouple

Ett histogram med 5000 slumpmässiga värden samplade från en skev gammafördelning ovan, och motsvarande histogram för medcouple-kärnvärdena nedan. Det faktiska medcouple är medianen för bottenfördelningen, markerad till 0,188994 med en gul linje.

Inom statistik är medcouple en robust statistik som mäter skevheten i en univariat fördelning . Den definieras som en skalad medianskillnad mellan vänster och höger halva av en distribution. Dess robusthet gör den lämplig för att identifiera extremvärden i justerade boxplots . Vanliga boxplotter klarar sig inte bra med snedfördelningar, eftersom de betecknar de längre osymmetriska svansarna som extremvärden. Med hjälp av medcouple kan morrhåren på en boxplot justeras för skevfördelningar och därmed få en mer exakt identifiering av extremvärden för icke-symmetriska fördelningar.

Som en sorts ordningsstatistik tillhör medparet klassen ofullständig generaliserad L-statistik . Liksom den vanliga medianen eller medelvärdet är medcouplet en icke-parametrisk statistik , så den kan beräknas för vilken fördelning som helst.

Definition

Följande beskrivning använder nollbaserad indexering för att harmonisera med indexeringen i många programmeringsspråk.

Låt $X:=\{x_{0}\geq x_{1}\geq \ldots \geq x_{n-1}\}$ vara ett ordnat urval av storlek $n$ , och låt $x_{m}$ vara medianen av $X$ . Definiera uppsättningarna

${\displaystyle X^{+}:=\{x_{i}~|~x_{i}\geq x_{m}\}} ,$ X

$X^{-}:=\{x_{j}~|~x_{j}\leq x_{m}\}$ ,

av storlekar $p:=|X^{+}|$ och $q:=|X^{-}|$ respektive. För $x_{i}^{+}\in X^{+}$ och $x_{j}^{-}\in X^{- }$ , definierar vi kärnfunktionen

$h(x_{i}^{+},x_{j}^{-}):={\begin{cases} \displaystyle {\frac {(x_{i}^{+}-x_{m})-(x_{m}-x_{j}^{-})}{x_{i}^{+}-x_{ j}^{-}}}&{\text{ if }}x_{i}^{+}>x_{j}^{-},\\\operatörsnamn {signum} (p-1-ij)&{ \text{ if }}x_{i}^{+}=x_{m}=x_{j}^{-},\end{cases}}$

där $\operatorname {signum}$ är teckenfunktionen .

Medparet är då median för uppsättningen

$\{h(x_{i}^{+},x_{j}^{-})~|~x_{i}^{+} \i X^{+}{\text{ och }}x_{j}^{-}\i X^{-}\}$ .

Med andra ord delar vi upp fördelningen i alla värden större eller lika med medianen och alla värden mindre än eller lika med medianen. Vi definierar en kärnfunktion vars första variabel är över de $p$ större värdena och vars andra variabel är över de $q$ mindre värdena. För det speciella fallet med värden kopplade till medianen, definierar vi kärnan med signum-funktionen . Medparet är då medianen över alla $pq$ -värden av $h(x_{i}^{+},x_{j}^{-})$ .

Eftersom medcouplet inte är en median som tillämpas på alla $(x_{i},x_{j})$ par, utan bara på de för vilka ${\displaystyle x_{i}^{+}\geq x_{m}\geq x_{j}^{-}} , den$ tillhör klassen ofullständig generaliserad L-statistik .

Egenskaper hos medparet

Medparet har ett antal önskvärda egenskaper. Några av dem ärvs direkt från kärnfunktionen.

Medcouple kärnan

Vi gör följande observationer om kärnfunktionen $h(x_{i}^{+},x_{j}^{-})$ :

1. Kärnfunktionen är platsinvariant. Om vi ​​adderar eller subtraherar något värde till varje element i exemplet $X$ ändras inte motsvarande värden för kärnfunktionen.
2. Kärnfunktionen är skalinvariant. Liknande skalning av alla element i exempel $X$ ändrar inte värdena för kärnfunktionen.

Dessa egenskaper ärvs i sin tur av medparet. Således är medcouple oberoende av medelvärdet och standardavvikelsen för en fördelning, en önskvärd egenskap för att mäta skevhet . För att underlätta beräkningen gör dessa egenskaper det möjligt för oss att definiera de två uppsättningarna

$Z^{+}:=\left.\left\{{\frac {x_{i}^{+}-x_{m}}{r}}~\right| ~x_{i}^{+}\in X^{+}\right\}$

$Z^{-}:=\left.\left\{{\frac {x_{j}^{-}-x_{m}}{r}}~\right| ~x_{j}^{-}\in X^{-}\right\}$

där $r=2\max _{0\leq i\leq n-1}|x_{i}|$ . Detta gör att mängden $Z:=Z^{+}\cup Z^{-}$ har ett intervall på högst 1, median 0, och behåller samma medcouple som ${\ displaystil X}$ .

För $Z$ reduceras medcouple-kärnan till

$h(z_{i}^{+},z_{j}^{-}):={\begin{cases}\displaystyle {\frac {z_{i}^{+} +z_{j}^{-}}{z_{i}^{+}-z_{j}^{-}}}&{\text{ if }}z_{i}^{+}>z_{j }^{-}\\\operatörsnamn {signum} (p-1-ij)&{\text{ if }}z_{i}^{+}=0=z_{j}^{-}\end{cases }}$

Genom att använda den nyligen ändrade och omskalade uppsättningen $Z$ kan vi observera följande.

3. Kärnfunktionen är mellan -1 och 1, det vill säga $|h(z_{i}^{+},z_{j}^{-})|\leq 1$ . Detta följer av den omvända triangelolikheten $|a|-|b|\leq |ab|$ med $a=z_{i}^{+}$ och $b=z_{j }^{-}$ och det faktum att $z_{i}^{+}\geq 0\geq z_{j}^{-}$ .
4. Medcouple kärnan $h(z_{i}^{+},z_{j}^{-})$ är icke-minskande i varje variabel. Detta kan verifieras av partiella derivator ${\frac {\partial h}{\partial z_{i}^{+}}}$ och ${ \frac {\partial h}{\partial z_{j}^{-}}}$ , båda icke-negativa, eftersom $z_{i}^{+}\geq 0\geq z_{j}^{-}$ .

Med egenskaperna 1, 2 och 4 kan vi alltså definiera följande matris ,

$H:=(h_{ij})=(h(z_{i}^{+},z_{j}^{-}))={\begin{pmatrix}h(z_{0}^{ +},z_{0}^{-})&\cdots &h(z_{0}^{+},z_{q-1}^{-})\\\vdots &\ddots &\vdots \\h (z_{p-1}^{+},z_{0}^{-})&\cdots &h(z_{p-1}^{+},z_{q-1}^{-})\end {pmatrix}}.$

Om vi ​​sorterar mängderna $Z^{+}$ och $Z^{-}$ i fallande ordning, så har matrisen $H$ sorterade rader och sorterade kolumner,

$H={\begin{pmatrix}h(z_{0}^{+},z_{0}^{-})&\geq &\cdots &\geq &h(z_{0}^{+} ,z_{q-1}^{-})\\\geq &&&&\geq \\\vdots &&\ddots &&\vdots \\\geq &&&&\geq \\h(z_{p-1}^{+} ,z_{0}^{-})&\geq &\cdots &\geq &h(z_{p-1}^{+},z_{q-1}^{-})\end{pmatrix}}.$

Medcouplet är då medianen för denna matris med sorterade rader och sorterade kolumner. Det faktum att raderna och kolumnerna är sorterade tillåter implementeringen av en snabb algoritm för att beräkna medcouple.

Robusthet

Nedbrytningspunkten är antalet värden som en statistik kan motstå innan den blir meningslös, dvs antalet godtyckligt stora extremvärden som datamängden $X$ kan ha innan värdet på statistiken påverkas. För medparet är nedbrytningspunkten 25 %, eftersom det är en median som tagits över paren $(x_{i},x_{j})$ så att $x_{i}\geq x_{m}\geq x_{j}$ .

Värderingar

Liksom alla mått på skevhet är medcouplet positivt för fördelningar som är sneda åt höger, negativt för fördelningar som är sneda åt vänster och noll för symmetriska fördelningar. Dessutom är värdena för medcouplet avgränsade av 1 i absolut värde.

Algoritmer för att beräkna medcouple

Innan vi presenterar medcouple-algoritmer minns vi att det finns ${\displaystyle O(n)}-$ algoritmer för att hitta medianen . Eftersom medcouple är en median, är vanliga algoritmer för median-finnande viktiga.

Naiv algoritm

Den naiva algoritmen för att beräkna medcouple är långsam. Det går vidare i två steg. Först konstruerar den medcouple-matrisen $H$ som innehåller alla möjliga värden för medcouple-kärnan. I det andra steget hittar den medianen för denna matris. Eftersom det finns $pq\approx {\frac {n^{2}}{4}}$ poster i matrisen i fallet när alla element i datamängden $X$ är unika, är den naiva algoritmens algoritmiska komplexitet $O(n^{2})$ .

Mer konkret fortsätter den naiva algoritmen enligt följande. Kom ihåg att vi använder nollbaserad indexering .

    
    
    
    
     
    
    
     function  naivt_medcouple(  vektor  X):  // X är en vektor med storlek n.  // Sortering i fallande ordning kan göras på plats i O(n log n) time  sort_decreasing  (X) xm := median(X) xscale := 2 * max(abs(X))  // Definiera övre och nedre centrerade och omskalade vektorer  // de ärver X:s egen minskande sortering  Zplus := [(x - xm)/xscale | x   i  X  så att  x >= xm] Zminus := [(x - xm)/xscale | x   i  X  så att  x <= xm] p := storlek(Zplus) q := storlek(Zminus)  // Definiera kärnfunktionen  som stänger  över Zplus och Zminus  funktion  h(i, j): a := Zplus[i] b := Zminus[j]  om  a == b:  returnera  signum  (p - 1 - i - j)  annars  :  returnera  (a + b) / (a ​​- b)  endif  endfunction  // O(n^2) operationer nödvändiga för att bilda denna vektor  H := [h(i, j) | i   i  [0, 1, ..., p - 1]  och  j  i  [0, 1, ..., q - 1]]  returnerar  median(H)  slutfunktion 

Det sista anropet till median på en vektor med storleken $O(n^{2})$ kan göras själv i $O(n^{2})$ operationer , därför är hela den naiva medcouple-algoritmen av samma komplexitet.

Snabb algoritm

Den snabba algoritmen överträffar den naiva algoritmen genom att utnyttja den sorterade naturen hos medcouple-matrisen $H$ . Istället för att beräkna alla poster i matrisen använder den snabba algoritmen K: ^te paralgoritmen från Johnson & Mizoguchi.

Det första steget av den snabba algoritmen fortsätter som den naiva algoritmen. Vi beräknar först de nödvändiga ingredienserna för kärnmatrisen, ${\displaystyle H=(h_{ij})} ,$ med sorterade rader och sorterade kolumner i fallande ordning. Istället för att beräkna alla värden för $h_{ij}$ , utnyttjar vi istället monotoniteten i rader och kolumner, via följande observationer.

Jämför ett värde med kärnmatrisen

Först noterar vi att vi kan jämföra vilken ${\displaystyle u} som helst$ med alla värden $h_{ij}$ av $H$ i $O(n)$ tid . Till exempel, för att bestämma alla $i$ och $j$ så att $h_{ij}>u$ har vi följande funktion:

             
         
         
         
         
           
         
         
           0
         
         
                    0
                  
             
                     
                     
             
             
             
                 
         
         
          
      funktion  större_h  (  kärna  h  ,  int  p  ,  int  q  ,  verklig  u  )  :  // h är kärnfunktionen, h(i,j) ger ith, j:te inmatningen av H  // p och q är antalet rader och kolumner av kärnmatrisen H  // vektor med storlek p  P  :  =  vektor  (  p  )  // indexering från noll  j  :  =  // med början från botten, beräkna [[supremum|minst övre gräns]] för varje rad  för  i  : =  p  -  1  ,  p  -  2  ,  ...,  1  ,  :  // sök på den här raden tills vi hittar ett värde mindre än u  medan  j  <  q  och  h  (  i  ,  j  )  >  u  :  j  :  =  j  +  1  samtidigt  // posten som föregår den vi just hittade är större än u  P  [  i  ]  :=  j  -  1  slutför  retur  P  slutfunktion 

Denna greater_h -funktion korsar kärnmatrisen från det nedre vänstra till det övre högra hörnet och returnerar en vektor $P$ av index som anger för varje rad var gränsen ligger mellan värden större än $u$ och de mindre än eller lika med $u$ . Denna metod fungerar på grund av den rad-kolumn sorterade egenskapen $H=(h_{ij})$ . Eftersom greater_h som mest beräknar $p+q$ värden för $h_{ij}$ , är dess komplexitet $O(n)$ .

Konceptuellt kan den resulterande $P$ -vektorn visualiseras som att upprätta en gräns på matrisen enligt följande diagram, där de röda posterna alla är större än $u$ :

Den symmetriska algoritmen för att beräkna värdena för $h_{ij}$ mindre än $u$ är mycket lika. Den fortsätter istället längs $H$ i motsatt riktning, från det övre högra hörnet till det nedre vänstra:

             
     
         
           
         
         
             
         
         
            0        
         
             
                0     
                     
             
             
             
                 
         
         
          
      funktion  less_h  (  kärna  h  ,  int  p  ,  int  q  ,  reell  u  )  :  // vektor med storlek p  Q  :  =  vektor  (  p  )  // sista möjliga radindex  j  :  =  q  -  1  // med början från toppen, beräkna [[infimum|största nedre gräns]] för varje rad  för  i  :=  ,  1  ,  ...,  p  -  2  ,  p  -  1  :  // sök på den här raden tills vi hittar ett värde som är större än u  medan  j  >=  och  h  (  i  ,  j  )  <  u  :  j  :  =  j  -  1  endwhile  // posten efter den vi nyss hittade är mindre än u  Q  [  i  ]  :=  j  +  1  endfor  return  Q  endfunction 

Denna nedre gräns kan visualiseras så, där de blå posterna är mindre än $u$ :

För varje $i$ har vi att ${\displaystyle P_{i}\geq Q_{i}} ,$ med strikt olikhet som endast förekommer för de rader som har värden lika med $u$ .

Vi har också att summorna

$\sum _{i=0}^{p-1}(P_{i}+1)~\qquad ~\ summa _{i=0}^{p-1}Q_{i}$

ange antalet element i $H$ som är större än $u$ , respektive antalet element som är större än eller lika med $u$ . Således ger denna metod också rangen för $u$ inom elementen $h_{ij}$ i $H$ .

Viktad median av radmedian

Den andra observationen är att vi kan använda den sorterade matrisstrukturen för att omedelbart jämföra vilket element som helst med minst hälften av posterna i matrisen. Till exempel är medianen för radmedianerna över hela matrisen mindre än den övre vänstra kvadranten i rött, men större än den nedre högra kvadranten i blått:

Mer generellt, med hjälp av gränserna som ges av vektorerna $P$ och $Q$ från föregående avsnitt, kan vi anta att vi efter några iterationer har bestämt positionen för medparet att ligga mellan den röda vänstra gräns och den blå högra gränsen:

De gula posterna anger medianen för varje rad. Om vi ​​mentalt omarrangerar raderna så att medianerna anpassas och ignorerar de kasserade posterna utanför gränserna,

vi kan välja en viktad median för dessa medianer, varje post viktad med antalet återstående poster på denna rad. Detta säkerställer att vi kan kassera minst 1/4 av alla återstående värden oavsett om vi måste kassera de större värdena i rött eller de mindre värdena i blått:

Varje radmedian kan beräknas i ${\displaystyle O(1)}-$ tid, eftersom raderna är sorterade, och den viktade medianen kan beräknas i ${\displaystyle O(n)}-$ tid, med hjälp av en binär sökning.

K: ^te parets algoritm

En visualisering av den snabba medcouple-algoritmen. Det börjar med en matris med sorterade rader och sorterade kolumner, där mörkare rutor är mindre än ljusare rutor. Vid varje iteration plockas den viktade medianen av radmedianerna i gult. Den jämförs sedan med resten av matrisen för att framställa kandidatröda övre och blå undre gränser. Algoritmen väljer sedan gränsen som är känd för att exkludera den globala matrismedianen, genom att beakta antalet poster som exkluderas av denna gräns (vilket är ekvivalent med att betrakta rangen för den gula posten). Algoritmen fortsätter sedan tills den gula viktade medianen av radmedianerna är exakt medparet, eller antalet kandidatposter är tillräckligt litet för att utföra en urvalssortering bland de återstående posterna.

Om man sätter samman dessa två observationer, fortskrider den snabba medcouple-algoritmen i stort sett enligt följande.

1. Beräkna de nödvändiga ingredienserna för medcouple-kärnfunktionen $h(i,j)$ med $p$ sorterade rader och $q$ sorterade kolumner.
2. Vid varje iteration, approximera medcouplet med den viktade medianen för radmedianerna.
3. Jämför denna preliminära gissning med hela matrisen för att erhålla höger och vänster gränsvektor $P$ respektive $Q$ . Summan av dessa vektorer ger oss också rangen för detta preliminära medpar.
   1. Om rangen för det preliminära medparet är exakt ${\displaystyle pq/2},$ sluta då. Vi har hittat medparet.
   2. I annat fall, kassera poster som är större än eller mindre än den preliminära gissningen genom att välja antingen $P$ eller $Q$ som den nya högra eller vänstra gränsen, beroende på vilken sida elementet i rang ${ \displaystyle pq/2}$ är in. Detta steg kasserar alltid minst 1/4 av alla återstående poster.
4. När antalet kandidatmedpar mellan höger och vänster gräns är mindre än eller lika med ${\displaystyle p} ,$ utför ett rangval bland de återstående posterna, så att rangordningen inom denna mindre kandidatuppsättning motsvarar $pq/2$ rang av medparet inom hela matrisen.

Den initiala sorteringen för att bilda ${\displaystyle h(i,j)}-$ funktionen tar $O(n\log n)$ tid. Vid varje iteration tar den viktade medianen $O(n)$ tid, såväl som beräkningarna av de nya preliminära $P$ och $Q$ vänster och höger gränser. Eftersom varje iteration kasserar minst 1/4 av alla återstående poster, kommer det att finnas som mest $O(\log n)$ iterationer. Alltså tar hela den snabba algoritmen $O(n\log n)$ tid.

Låt oss återge den snabba algoritmen mer i detalj.

    
    
     
    
    
    
    
    
    
    
    
         
        
        
   
      function  medcouple(  vektor  X):  // X är en vektor med storlek n  // Beräkna initiala ingredienser som för  det naiva medcouple  sort_decreasing  (X) xm := median(X) xscale := 2 * max(abs(X)) Zplus := [(x - xm)/xskala | x   i  X  så att  x >= xm] Zminus := [(x - xm)/xscale | x   i  X  så att  x <= xm] p := storlek(Zplus) q := storlek(Zminus)  funktion  h(i, j): a := Zplus[i] b := Zminus[j]  om  a == b:  return  signum  (p - 1 - i - j)  else  :  return  (a + b) / (a ​​- b)  endif  endfunction  // Begin Kth par algoritm (Johnson & Mizoguchi)  // De initiala vänster- och högergränserna, två vektorer med storleken p  L := [0, 0, ..., 0] R := [q - 1, q - 1, ..., q - 1]  // antal poster till vänster om den vänstra gränsen  Ltotal := 0  // antal poster till vänster om den högra gränsen  Rtotal := p*q  // Eftersom vi indexerar från noll är medcouple index en  // mindre än dess rang.  medcouple_index :=  floor  (Rtotal / 2)  // Iterera medan antalet poster mellan gränserna är  // större än antalet rader i matrisen.  medan  Rtotal - Ltotal > p:  // Beräkna radmedianer och deras associerade vikter, men hoppa över  // alla rader som redan är tomma.  middle_idx := [i | i   i  [0, 1, ..., p - 1]  så  att  L[i] <= R[i]] rad_medianer := [h(i,  golv  ((L[i] + R[i])/ 2) | i  i  middle_idx] vikter := [R[i] - L[i] + 1 | i  i  middle_idx] WM :=  viktad median  (rad_medianer, vikter)  // Nya tentativa höger- och vänstergränser  P :=  greater_h  ( h, p, q, WM) Q :=  mindre_h  (h, p, q, WM) Ptotal := summa(P) + storlek(P) Qtotal := summa(Q) // Bestäm vilka  poster som ska kasseras, eller om vi har hittat medcouple  om  medcouple_index <= Ptotal - 1: R := P Rtotal := Ptotal  else  :  if  medcouple_index > Qtotal - 1: L := Q Ltotal := Qtotal  else  : // Hittade medcouple, rank of the viktad median är lika med medcouple index  return  WM  endif  endif  endwhile  // Hittade inte medcouple, men det finns väldigt få preliminära poster kvar := [h(i, j) | i  in  [0, 1, ..., p - 1], j  i  [L[i], L[i] + 1, ..., R[i]]  så  att  L[i] <= R[i] ]  // Välj medparet efter rangordning bland de återstående poster  medcouple :=  select_nth  (resterande, medcouple_index - Ltotal)  returnera  medcouple  endfunction 

I verklig användning måste algoritmen också ta hänsyn till fel som uppstår från flyttalsaritmetik med ändlig precision . Till exempel bör jämförelserna för medcouple kernel-funktionen göras inom maskin epsilon , liksom ordningsjämförelserna i greater_h och less_h -funktionerna .

Programvara/källkod

• Den snabba medcouple-algoritmen är implementerad i R :s robustbase-paket.
• Den snabba medcouple-algoritmen är implementerad i en C-förlängning för Python i Robustats Python-paketet .
• En GPL-ad C++ - implementering av den snabba algoritmen , härledd från R-implementeringen.
• En Stata- implementering av den snabba algoritmen .
• En implementering av den naiva algoritmen i Matlab (och därmed GNU Octave ).
• Den naiva algoritmen är också implementerad för Python- paketets statistikmodeller .

Se även

• Robust statistik
• Skevhet
• Justerade boxplotter