## General functions
rpareto <- function(n, m, a) (m / (1-runif(n))**a)
MAD <- function(x) median(abs(x - median(x)))
M <- function(x) (0.6745 * (x - median(x))/MAD(x))

## Try on an example
N <- 200
samples <- rpareto(N, 1, 1)
D <- abs(M(samples))
cutoff <- 3.5
print(paste("Fraction of values marked as outliers:", sum(D > cutoff)/N))
print(paste("Lowest outlier:", min(samples[D > cutoff])))
E <- sum(samples)
Ee <- sum(samples[D > cutoff])
print(paste("Contributions of outliers to expectation:", Ee/E))