Это продолжение, но также другой вопрос моего предыдущего .
Я читал в Википедии, что « Средне-несмещенный оценщик минимизирует риск по отношению к функции потери абсолютного отклонения, как это наблюдал Лаплас ». Тем не менее, мои результаты моделирования Монте-Карло не поддерживают этот аргумент.
Я предполагаю выборку из логарифмически нормального населения, , где, и \ sigma - это среднее значение логарифма и log-sd, \ beta = \ exp (\ mu) = 50μ σ β = exp ( μ ) = 50
Среднегеометрическая оценка - это несмещенная по медиане оценка для медианы населения ,
где и - среднее значение log, а log-sd, и - MLE для и .
В то время как исправленная средняя геометрическая оценка является средней несмещенной оценкой для медианы населения.
Я генерирую образцы размера 5 повторно из LN . Номер репликации составляет 10000. Среднее абсолютное отклонение, которое я получил, составляет 25,14 для среднего геометрического и 22,92 для исправленного среднего геометрического. Почему?
Кстати, оцененные средние абсолютные отклонения составляют 18,18 для среднего геометрического и 18,58 для скорректированного среднего геометрического.
Я использовал скрипт R здесь:
#```{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
#```
set.seed
. 3.) Не всегда доверяйте Википедии - обратите внимание, чем ваш цитируемый текст (из статьи «Медиана») отличается от этой другой статьи 4. Википедия . Ваш код R - полный беспорядок - посмотрите руководство по стилю Google R для некоторых хороший стиль руководства.Ответы:
Если мы выберем оценщик по критерию, который минимизирует ожидаемую абсолютную ошибку от истинного значения αα+ α
мы требуем
что эквивалентно . Таким образом, является медианой, следующей за Лапласом в 1774 году.п( α > α+) = 1 / 2 α+
Если у вас возникли проблемы с R, задайте их в другом вопросе о переполнении стека.
источник