R-Code für den verallgemeinerten (Extreme Studentized Deviate) ESD-Outlier-Test

Question

Hier ist der Test, an dem ich interessiert bin: http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm

Wie kann ich diesen Code in eine Funktion umwandeln, die einen Vektor numerischer Werte akzeptiert und einen logischen Vektor zurückgibt, der angibt, welche Datenpunkte entfernt werden sollen?

Ich habe unten versucht, dies zu tun, aber ich komme nicht weiter, weil sich die sortierte Vektor zum Zurückgeben nicht mit den Eingabedaten des Vektors deckt.

# Eingabedaten
y = c(-0.25, 0.68, 0.94, 1.15, 1.20, 1.26, 1.26,
      1.34, 1.38, 1.43, 1.49, 1.49, 1.55, 1.56,
      1.58, 1.65, 1.69, 1.70, 1.76, 1.77, 1.81,
      1.91, 1.94, 1.96, 1.99, 2.06, 2.09, 2.10,
      2.14, 2.15, 2.23, 2.24, 2.26, 2.35, 2.37,
      2.40, 2.47, 2.54, 2.62, 2.64, 2.90, 2.92,
      2.92, 2.93, 3.21, 3.26, 3.30, 3.59, 3.68,
      4.30, 4.64, 5.34, 5.42, 6.01)

## Normalverteilungsdiagramm generieren.
qqnorm(y)

removeoutliers = function(dfinputcol) {

  y = as.vector(dfinputcol)

  ## Funktion zum Berechnen der Teststatistik erstellen.
  rval = function(y){
    ares = abs(y - mean(y))/sd(y)
    df = data.frame(y, ares)
    r = max(df$ares)
    list(r, df)}

  ## Werte und Vektoren definieren.
  n = length(y)
  alpha = 0.05
  lam = c(1:10)
  R = c(1:10)

  ## Teststatistik berechnen, bis r=10 Werte aus der Stichprobe entfernt wurden.
  for (i in 1:10){

    if(i==1){
      rt = rval(y)
      R[i] = unlist(rt[1])
      df = data.frame(rt[2])
      newdf = df[df$ares!=max(df$ares),]}

    else if(i!=1){
      rt = rval(newdf$y)
      R[i] = unlist(rt[1])
      df = data.frame(rt[2])
      newdf = df[df$ares!=max(df$ares),]}

    ## Kritischen Wert berechnen.
    p = 1 - alpha/(2*(n-i+1))
    t = qt(p,(n-i-1))
    lam[i] = t*(n-i) / sqrt((n-i-1+t**2)*(n-i+1))

  }
  ## Ergebnisse ausgeben.
  newdf = data.frame(c(1:10),R,lam)
  names(newdf)=c("Ausreißer","Teststatistik", "Kritischer Wert")

  # Festlegen, wie viele Ausreißer entfernt werden sollen
  toremove = max(newdf$Ausreißer[newdf$Teststatistik > newdf$Kritischer Wert])

  # Vektor gleicher Größe wie Eingabevektor erstellen
  logischerVektorWennsollentferntwerden = logical(length=length(y))

  # aber wie bestimmt man, welche Ausreißer entfernt werden sollen?
  # die größten Datenpunkte als Ausreißer zum Entfernen festlegen... können jedoch in einigen Datensätzen die kleinsten sein..
  logischerVektorWennsollentferntwerden = replace(logischerVektorWennsollentferntwerden, tail(sort(y), toremove), TRUE)

  return (logischerVektorWennsollentferntwerden)
}

# dies sollte 3 Datenpunkte auf TRUE gesetzt haben.. hat aber nur 2 und das sind nicht die richtigen
output = removeoutliers(y)
length(output[output==T])

Answer 1

Sie können winsorize in der Bibliothek robustHD verwenden

library('robustHD')

set.seed(1234) 
x <- rnorm(10) 
x[1] <- x[1] * 10 
x[2] <- x[2] * 11
x[10] <- x[10] * 10
x
[1] -12.0706575   3.0517217   1.0844412  -2.3456977   0.4291247   0.5060559  -0.5747400  -0.5466319  -0.5644520  -8.9003783
boxplot(x)

y <- winsorize(x)
y
 [1] -4.5609058  3.0517217  1.0844412 -2.3456977  0.4291247  0.5060559 -0.5747400 -0.5466319 -0.5644520 -4.5609058
boxplot(y)

Wenn Sie ein DataFrame oder einen Vektor haben, können Sie sapply verwenden, um die Winsorize-Funktion auszuführen. Weitere Informationen zu dieser Bibliothek finden Sie unter diesem Link http://cran.r-project.org/web/packages/robustHD/index.html

Answer 2

Ich denke, es steht auf der Seite, die du gegeben hast (nicht genau, aber in zwei Sätzen):

Entferne die r Beobachtungen, die |x_i - Mittelwert(x)| maximieren

Du erhältst die Daten ohne Ausreißer einfach, indem du die r Werte entfernst, die über der Differenz liegen, mit:

y[abs(y-mean(y)) < sort(abs(y-mean(y)),decreasing=TRUE)[toremove]]

Du benötigst nicht die letzten beiden Zeilen deines Codes. Übrigens, du kannst direkt berechnen:

toremove = which(newdf$TestStat > newdf$CriticalVal)

Um es etwas zu vereinfachen, könnte die endgültige Funktion folgendermaßen aussehen:

# Berechne den kritischen Wert für den ESD-Test
esd.critical <- function(alpha, n, i) {
    p = 1 - alpha/(2*(n-i+1))
    t = qt(p,(n-i-1))
    return(t*(n-i) / sqrt((n-i-1+t**2)*(n-i+1)))
}

removeoutliers = function(y) {

    ## Definition der Werte und Vektoren.
    y2 = y
    n = length(y)
    alpha = 0.05
    toremove = 0

    ## Berechnung der Teststatistik, bis r=10 Werte aus der Stichprobe entfernt wurden.
    for (i in 1:10){
        if(sd(y2)==0) break
        ares = abs(y2 - mean(y2))/sd(y2)
        Ri = max(ares)
        y2 = y2[ares!=Ri]

        ## Berechne kritischen Wert.
        if(Ri>esd.critical(alpha,n,i))
            toremove = i
    }

    # Werte, die behalten werden sollen
    if(toremove>0)
        y = y[abs(y-mean(y)) < sort(abs(y-mean(y)),decreasing=TRUE)[toremove]]

    return (y)
}

was folgendes zurückgibt:

> removeoutliers(y)
 [1] -0.25  0.68  0.94  1.15  1.20  1.26  1.26  1.34  1.38  1.43  1.49
[12]  1.49  1.55  1.56  1.58  1.65  1.69  1.70  1.76  1.77  1.81  1.91
[23]  1.94  1.96  1.99  2.06  2.09  2.10  2.14  2.15  2.23  2.24  2.26
[34]  2.35  2.37  2.40  2.47  2.54  2.62  2.64  2.90  2.92  2.92  2.93
[45]  3.21  3.26  3.30  3.59  3.68  4.30  4.64