您还可以使用重尾 Lambert W x F 分布估计柯西的位置（免责声明：我是作者。）因为两者都是对称的$c$（柯西的位置）和$\mu_x$（输入 X ~ F 的平均值），分别。事实上，我在论文中给出了一个估计 Cauchy 位置的示例，并比较了user777 建议的累积样本平均估计值。

对于 F 是正态分布和$\alpha = 1$, 变换后的随机变量$Y = func(X, \delta)$减少到 Tukey 的 h 分布。为了$\delta = 0$它们是正态分布；为了$\delta > 0$他们有更重的尾巴。Lambert W x F 分布的优点是您也可以再次从非正态返回到正态；即，您可以估计参数和Gaussianize()数据。

在 R 中，您可以使用LambertW 包模拟、估计、绘制等多个 Lambert W x F 分布。

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

现在让我们看看模拟

# 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

正如我们之前知道的那样，这mean是一个糟糕的估计器，所以我们将从图中删除它。

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

它们看起来都非常接近（median并且IGMM稍微更糟）。

柯西分布具有无限均值和无限方差。由于这个事实，大数定律和中心极限定理不适用。

此演示旨在让您直观地了解当这些样本是 Cauchy 时向样本添加额外观察时会发生什么。最终，您从分布中得出一个值，该值相对于其他值非常大，以至于它“消除”了恢复均值的影响。

增加样本量不会使均值“趋向于”柯西分布的真实位置。作为演示，编写一个程序来计算一个大数$n$的柯西偏差。第一个样本的平均值$1...i$英石$i\le n$偏差将在非常小的值和非常大的值之间剧烈波动。您可以在这些运行均值与用于计算运行均值的偏差数量的图中轻松看到这一点。

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

现在再添加 1000 个观察值x，看看会发生什么。

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

运行均值似乎仍未恢复$0$很快......似乎当它接近时，它会产生一个非常非常大的偏差，所以平均值会“跳”离分布的位置。