Remedian

Remedian

The function used here is a direct translation of the FORTRAN code in the following article. Rousseeuw, P. J., & Bassett Jr, G. W. (1990). The remedian: A robust averaging method for large data sets. Journal of the American Statistical Association, 85(409), 97-104.

remedian <- function(b = 15, k = 4, x){
  a1 <- c()
  a2 <- c()
  a3 <- c()
  a4 <- c()
  ll = 0
  for(i in 1:b){
    for(j in 1:b){
      for(k in 1:b){
        for(l in 1:b){
          ll <- ll + 1
          a1[l] <- x[ll]
        }
        a2[k] <- median(a1)
      }
      a3[j] <- median(a2)
    }
    a4[i] <- median(a3)
  }
  
  return(median(a4))
}

Suppose x is a vector of n observations of a standard normal random variable, where n = 50625, b = 15, k = 4, \(n = b^k\). Calculate remedian 100 times.

rem <- c()
me <- c()
for(h in 1:100){
x <- rnorm(50625)
rem[h] <- remedian(15,4,x)
me[h] <- median(x)
}

A summary of the results.

summary(rem)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.026203 -0.003676  0.001649  0.001693  0.007364  0.035430

plot(rem~me)

summary(lm(rem~me))

## 
## Call:
## lm(formula = rem ~ me)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0243817 -0.0055634 -0.0003713  0.0056328  0.0279134 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0013654  0.0008716   1.566     0.12    
## me          0.9528773  0.1676976   5.682 1.37e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.008697 on 98 degrees of freedom
## Multiple R-squared:  0.2478, Adjusted R-squared:  0.2401 
## F-statistic: 32.29 on 1 and 98 DF,  p-value: 1.367e-07