Comparing estimators of location of the Cauchy distribution

Question

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

Answer 1

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.

I also suggest reading Why does the Cauchy distribution have no mean?

Answer 2

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