## The proportion of outliers.
alpha = 0.05

## The variance of the outlier component.
tau2 = 100

## The variance of the non-outlier component.
sig2 = 1

## Number of replications.
nrep = 1e3

## Sample size
n = 50

## Inidicators of where we have outliers.
A = array((runif(n*nrep) < alpha), c(nrep,n))

## The outliers.
B = array(sqrt(tau2)*rnorm(n*nrep), c(nrep,n))

## The non-outliers.
C = array(sqrt(sig2)*rnorm(n*nrep), c(nrep,n))

## These are the contaminated normal draws.
X = A*B + (1-A)*C

## Sort the rows.
X = apply(X, 1, sort)
X = t(X) ## R peculiarity

## Use trimmed means for a range of values of k.
MSE = NULL
for (k in (0:20))
{
  TX = X[,(k+1):(n-k)]
  TM = apply(TX, 1, mean)
  MSE[k] = mean(TM^2)
}

RMSE = sqrt(MSE)