Apakah estimator median-bias meminimalkan rata-rata penyimpangan absolut?

14

Ini adalah tindak lanjut tetapi juga pertanyaan yang berbeda dari pertanyaan saya sebelumnya .

Saya membaca di Wikipedia bahwa " Estimator yang tidak bias rata-rata meminimalkan risiko sehubungan dengan fungsi kerugian absolut-deviasi, seperti yang diamati oleh Laplace ." Namun, hasil simulasi Monte Carlo saya tidak mendukung argumen ini.

Saya menganggap sampel dari populasi log-normal, , di mana, dan adalah log-mean dan log-sd,μ σ β = exp ( μ ) = 50X1,X2,...,XNLN(μ,σ2)μσβ=exp(μ)=50

Penduga rata-rata geometrik adalah penduga rata-rata median untuk median populasi ,exp(μ)

β^GM=exp(μ^)=exp(log(Xi)N)LN(μ,σ2/N) mana, dan adalah log-mean dan log-sd, dan adalah MLE untuk dan .μσμ^σ^μσ

Sementara penduga rata-rata geometrik yang diperbaiki adalah penduga rata-rata yang tidak bias untuk median populasi.

β^CG=exp(μ^σ^2/2N)

Saya membuat sampel ukuran 5 berulang kali dari LN . Jumlah replikasi adalah 10.000. Penyimpangan absolut rata-rata yang saya dapatkan adalah 25,14 untuk penduga rata-rata geometrik dan 22,92 untuk rata-rata geometrik terkoreksi. Mengapa?(log(50),log(1+22))

BTW, estimasi penyimpangan absolut rata-rata adalah 18,18 untuk rata-rata geometrik dan 18,58 untuk penduga rata-rata geometrik yang dikoreksi.

Skrip R yang saya gunakan ada di sini:

#```{r stackexchange}
#' Calculate the geomean to estimate the lognormal median.
#'
#' This function Calculate the geomean to estimate the lognormal
#' median.
#'
#' @param x a vector.
require(plyr)
GM <- function(x){
    exp(mean(log(x)))
}
#' Calculate the bias corrected geomean to estimate the lognormal
#' median.
#'
#' This function Calculate the bias corrected geomean using the
#' variance of the log of the samples, i.e., $\hat\sigma^2=1/(n-1)
# \Sigma_i(\Log(X_i)-\hat\mu)^2$
#'
#' @param x a vector.
BCGM <- function(x){
y <- log(x)
exp(mean(y)-var(y)/(2*length(y)))
}
#' Calculate the bias corrected geomean to estimate the lognormal
#' median.
#'
#' This function Calculate the bias corrected geomean using
#' $\hat\sigma^2=1/(n)\Sigma_i(\Log(X_i)-\hat\mu)^2$
#'
#' @param x a vector.
CG <- function(x){
y <- log(x)
exp(mean(y)-var(y)/(2*length(y))*(length(y)-1)/length(y))
}

############################

simln <- function(n,mu,sigma,CI=FALSE)
{
    X <- rlnorm(n,mu,sigma)
    Y <- 1/X
    gm <- GM(X)
    cg <- CG(X)
    ##gmk <- log(2)/GM(log(2)*Y) #the same as GM(X)
    ##cgk <- log(2)/CG(log(2)*Y)
    cgk <- 1/CG(Y)
    sm <- median(X)
    if(CI==TRUE) ci <- calCI(X)
    ##bcgm <- BCGM(X)
    ##return(c(gm,cg,bcgm))
    if(CI==FALSE) return(c(GM=gm,CG=cg,CGK=cgk,SM=sm)) else return(c(GM=gm,CG=cg,CGK=cgk,CI=ci[3],SM=sm))
}
cv <-2
mcN <-10000
res <- sapply(1:mcN,function(i){simln(n=5,mu=log(50),sigma=sqrt(log(1+cv^2)), CI=FALSE)})
sumres.mad <- apply(res,1,function(x) mean(abs(x-50)))
sumres.medad <- apply(res,1,function(x) median(abs(x-50)))
sumres.mse <- apply(res,1,function(x) mean((x-50)^2))
#```

#```{r eval=FALSE}
#> sumres.mad
      GM       CG      CGK       SM 
#25.14202 22.91564 29.65724 31.49275 
#> sumres.mse
      GM       CG      CGK       SM 
#1368.209 1031.478 2051.540 2407.218 
#```
Zhenglei
sumber
1
1.) "10.000" terlalu kecil untuk pertanyaan Anda - coba "250.000" (atau lebih). 2.) Jika Anda menjalankan simulasi Monte Carlo dan mendapatkan hasil yang tampak aneh, coba ganti dengan seed set.seed. 3.) Jangan selalu percaya pada Wikipedia - perhatikan bagaimana teks yang Anda kutip (dari artikel "Median") berbeda dari artikel Wikipedia lainnya 4.) Kode R Anda berantakan total - lihat Panduan Gaya R Google untuk beberapa pedoman gaya yang baik.
Steve S

Jawaban:

4

α+α

E=<|α+α|>=α+(α+α)f(α)dα+α+(αα+)f(α)dα

kami membutuhkan

dEdα+=α+f(α)dαα+f(α)dα=0

P(α>α+)=1/2α+

Jika Anda mengalami masalah dengan R, silakan tanyakan dalam pertanyaan lain di Stack Overflow

Keith
sumber
Secara teoritis, saya pikir itu benar. Namun, saya bingung dengan hasil simulasi R yang tidak mendukung pernyataan ini seperti yang diharapkan.
Zhenglei
2
Saya seorang Ilmuwan Data / Fisika sehingga belum pernah melihat garis R. Seperti yang saya sarankan dalam pertanyaan, jika itu adalah masalah kode Anda harus menanyakannya di Stack Overflow dan Anda akan mendapatkan lebih banyak perhatian. Namun, jawaban di atas adalah benar kecuali jika Anda ingin menguraikan bagaimana generalisasi ke estimator median-tidak memihak. Untuk lebih jelasnya lihat halaman 172 dari buku ET Jaynes. Probability theory ISBN 978-0-521-59271-0.
Keith
Terima kasih banyak atas jawaban Anda. Ini bukan masalah pengkodean. Saya hanya ingin melakukan simulasi untuk menunjukkan bahwa estimator median-bias akan meminimalkan deviasi absolut yang diharapkan. Saya belum menerima jawabannya karena saya terutama bingung tentang langkah simulasi. Saya menerapkannya dalam R tetapi simulasi dapat dilakukan dalam Matlab atau Python atau bahasa lain.
Zhenglei
2
@Keith maaf untuk matematika saya yang lemah, tetapi bisakah Anda menunjukkan lebih detail tentang bagaimana Anda mendapatkan ekspektasi tersebut?
AdamO