## 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))