Cauchy分布のパラメタ推定

Cauchy分布の乱数を100個生成し，その中心の値をどうすれば正確に求められるかを試してみる。

sim = function() {
    r = rcauchy(100)
    m1 = mean(r)                  # 平均
    m2 = median(r)                # 中央値
    f = function(par) {
        -sum(dcauchy(r, par[1], par[2], log=TRUE))
    }
    m3 = optim(c(0,1), f)$par[1]  # 最尤推定値
    c(m1, m2, m3)
}

r = replicate(100000, sim())

それぞれの誤差を求めてみる：

sqrt(mean(r[1,]^2))
[1] 451.953
sqrt(mean(r[2,]^2))
[1] 0.158399
sqrt(mean(r[3,]^2))
[1] 0.1432506
max(abs(r[1,]))
[1] 131453.7
max(abs(r[2,]))
[1] 0.7711901
max(abs(r[3,]))
[1] 0.677196

最後の方法が一番正確だが，中央値とは僅差である。