# Comparing estimators of location of the Cauchy distribution

I'm comparing the following 4 estimators of location of the Cauchy distribution: Let $x_{1},..x_{n}$ be observations and $l$ be the log likelihood function.

$x=median(x_{1},..x_{n})$, $y=x+\frac{l'(x)}{n/2}$, $z=mean(x_{1},..x_{n})$ and $s=MLE$

I did simulations in R of the mean squared error and probability coverages. I found that the mean was a very bad estimator. The best estimator was $y$ in terms of giving the smallest MSE and the greatest probability coverage. $x$ and $z$ were pretty good estimators as well and in fact when I increased $n$ large enough $x$ and $s$ gave the same results. I know the reason $x$ is a bad estimator. What are the mathematical reasons for my other findings?

Thanks

• And how do you define the location of the distribution? Mar 19 '14 at 8:26
• Note the cauchy distribution have infinite mean too, so that is the reason why it performs poorly. Mar 19 '14 at 8:27
• Your question has a typo. You say the mean performs badly, then later you say that z was a pretty good estimator... but that's the mean. Mar 19 '14 at 8:29
• @mpiktas Yes I meant s, not z. Infinite variance implies infinite mean. So can you please help with the theoretical reasons for the results other than z? Thanks Mar 19 '14 at 10:42
• The mean of iid Cauchy variables also has a Cauchy distribution, whence it has infinite variance. This means that no matter how large the sample becomes, its mean will never become a better estimate of central tendency ("better" in a least-squares sense) than any single observation and it will have infinite quadratic loss. That's pretty bad.
– whuber
Mar 19 '14 at 14:57

Cauchy distributions have infinite mean and infinite variance. Because of this fact, the laws of large numbers and central limit theorems does not apply.

This demonstration is designed to give some intuition about what happens as you add additional observations to your sample when those samples are Cauchy. Eventually, you draw a value from the distribution which is so large relative to the other values that it "washes out" the effect of reverting to the mean.

Increased sample size will not make the mean "tend toward" the true location of the Cauchy distribution. For a demonstration, write a program to compute a large number $$n$$ of Cauchy deviates. The mean of the sample for the first $$1...i$$ s.t. $$i\le n$$ deviates will wildly oscillate between very small and very large values. You can see this easily in a plot of those running means versus number of deviates used to compute the running mean.

x   <- rcauchy(1000)
y   <- NULL
for(i in 1:length(x)){
y[i]    <- mean(x[1:i])
}
plot(1:length(x),y, type="l")


Now add another 1000 observations to x and see what happens.

x   <- c(x, rcauchy(1000))
y   <- NULL
for(i in 1:length(x)){
y[i]    <- mean(x[1:i])
}
plot(1:length(x),y, type="l")


The running mean still doesn't appear to be returning to $$0$$ very quickly... it almost seems as if when it gets close, it a very, very large deviate will be drawn, so the mean will "jump" away from the location of the distribution.

You can also estimate the location of the Cauchy using heavy-tail(s) Lambert W x F distributions (Disclaimer: I am the author.) since both are symmetric around $$c$$ (location of Cauchy) and $$\mu_x$$ (mean of the input X ~ F), respectively. In fact, I give an example of estimating the location of a Cauchy in the paper and compare the cumulative sample average estimates as suggested by user777.

For F being the Normal distribution and $$\alpha = 1$$, the transformed random variable $$Y = func(X, \delta)$$ reduces to Tukey's h distribution. For $$\delta = 0$$ they are the Normal distribution; for $$\delta > 0$$ they have heavier tails. The nice property of Lambert W x F distributions is that you can also go back from non-normal to Normal again; i.e., you can estimate parameters and Gaussianize() your data.

In R you can simulate, estimate, plot, etc. several Lambert W x F distributions with the LambertW package.

library(LambertW)
library(MASS)  # for fitdistr()

LogLikCauchy <- function(loc, x.sample) {
# sum(dcauchy(x.sample, location = loc, scale = 1, log = TRUE))
nn <- length(x.sample)
return(- nn * log(pi) - sum((log(1 + (x.sample - loc)^2))))
}

DerivLogLikCauchy <- function(loc, x.sample) {
return(sum(1 / (1 + (x.sample - loc)^2) * 2 * (x.sample - loc)))
}

LocationEstimators <- function(x.sample) {
nn <- length(x.sample)
out <-
c(mean = mean(x.sample),
median = median(x.sample),
mle.cauchy = suppressWarnings(fitdistr(x.sample,
"cauchy")$$est["location"])) out <- c(out, median.loglik = out["median"] + DerivLogLikCauchy(out["median"], x.sample) / (nn / 2)) # Lambert W x Gaussian estimates for heavy tails ('h') igmm.tau <- LambertW::IGMM(x.sample, "h")$$tau
beta.hat <- igmm.tau[1:2]
names(beta.hat) <- c("mu", "sigma")
mle.lambertw <- LambertW::MLE_LambertW(x.sample, distname = "normal",
theta.init = LambertW::tau2theta(igmm.tau,
beta = beta.hat),
type = "h",
return.estimate.only = TRUE)
out <- c(out, igmm.tau["mu_x"], mle.lambertw["mu"])
names(out)[3:6] <-
c("median.loglik", "mle.cauchy", "igmm.LambertW", "mle.LambertW")
return(out)
}


Now let's look at simulations

# simulate and look at bias, std dev, and MSE
nsim <- 1000
num.samples <- 100
set.seed(nsim)
est <- t(replicate(nsim,
LocationEstimators(rcauchy(num.samples))))

colMeans(est)
##          mean        median median.loglik    mle.cauchy igmm.LambertW
##      -0.43373       0.00255       0.00306       0.00237      -0.00326
##  mle.LambertW
##       0.00321

apply(est, 2, sd)
##          mean        median median.loglik    mle.cauchy igmm.LambertW
##        29.183         0.156         0.145         0.144         0.221
##  mle.LambertW
##         0.146

# RMSE (since true location = 0)
sqrt(colMeans(est^2))
##          mean        median median.loglik    mle.cauchy igmm.LambertW
##        29.171         0.156         0.145         0.144         0.221
##  mle.LambertW
##         0.146


As we knew beforehand, the mean is a bad estimator, so we'll remove it from the plots.

library(ggplot2)
library(reshape2)
theme_set(theme_bw(18))
est.m <- melt(est)
colnames(est.m) <- c("sim.id", "estimator", "value")

# remove 'mean' for good scaling in plots
est.m <- subset(est.m, estimator != 'mean')
ggplot(est.m,
aes(estimator, value, fill = estimator)) +
geom_violin() +
geom_hline(yintercept = 0, size = 1, linetype = "dashed",
colour = "blue") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 90))


They all seem pretty close to each other (with median and IGMM being slightly worse).