¿Un estimador imparcial medio minimiza la desviación absoluta media?

14

Este es un seguimiento, pero también una pregunta diferente de la anterior .

Leí en Wikipedia que " Un estimador imparcial mediano minimiza el riesgo con respecto a la función de pérdida de desviación absoluta, como lo observó Laplace ". Sin embargo, mis resultados de simulación de Monte Carlo no respaldan este argumento.

Asumo una muestra de una población-log normal, , donde, μ y σ son el log-mean y log-sd, β = exp ( μ ) = 50X1,X2,...,XNLN(μ,σ2)μσβ=exp(μ)=50

El estimador de la media geométrica es un estimador imparcial medio para la mediana de la población ,exp(μ)

, donde,μyσson la media logarítmica log-sd, μ y σ son los MLEs paraμyσ.β^GM=exp(μ^)=exp(log(Xi)N)LN(μ,σ2/N)μσμ^σ^μσ

Mientras que un estimador de media geométrica corregida es un estimador imparcial de media para la mediana de la población.

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

Genero muestras de tamaño 5 repetidamente desde el LN . El número de replicación es 10,000. Las desviaciones absolutas promedio que obtuve son 25.14 para el estimador de la media geométrica y 22.92 para la media geométrica corregida. ¿Por qué?(log(50),log(1+22))

Por cierto, las desviaciones absolutas medias estimadas son 18.18 para la media geométrica y 18.58 para el estimador de la media geométrica corregida.

El script R que utilicé está aquí:

#```{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
fuente
1
1.) "10,000" es demasiado pequeño para su pregunta - intente con "250,000" (o más). 2.) Si ejecuta una simulación de Monte Carlo y obtiene un resultado que parece extraño, intente cambiar la semilla con set.seed. 3.) No confíes siempre en Wikipedia: observa cómo tu texto citado (del artículo "Mediano") difiere de este otro artículo de Wikipedia 4.) Tu código R es un desastre total: consulta la Guía de estilo R de Google para obtener algunos buenas pautas de estilo.
Steve S

Respuestas:

4

α+α

E=<|α+α|>=α+(α+-α)F(α)reα+α+(α-α+)F(α)reα

necesitamos

remireα+=-α+F(α)reα-α+F(α)reα=0 0

que es equivalente a PAG(α>α+)=1/ /2. Entoncesα+ se muestra que es la mediana que sigue a Laplace en 1774.

Si tiene problemas con R, hágalo en otra pregunta sobre Stack Overflow

Keith
fuente
Teóricamente, creo que es correcto. Sin embargo, estoy confundido por los resultados de la simulación R que no respaldan esta declaración como se esperaba.
Zhenglei
2
Soy un científico de datos / físico, por lo que nunca he visto una línea de R. Como sugerí en la pregunta, si se trata de un problema de código, debe preguntarlo en Stack Overflow y obtendrá mucha más atención. Sin embargo, la respuesta anterior es correcta a menos que desee detallar cómo se generaliza a un estimador imparcial medio. Para más detalles ver la página 172 del libro ET Jaynes Teoría de la probabilidad ISBN 978-0-521-59271-0.
Keith
Thank you a lot for your answer. It is not a coding issue. I just want to do simulations to show that a median-unbiased estimator will minimize the expected absolute deviation. I haven't accepted the answer because I am mainly confused about the simulation step. I implemented it in R but simulations could be done in Matlab or Python or any other languages.
Zhenglei
2
I suspect the issue is that you are dealing with an approximation which works as N -> pero tienes 10,000 y 5 que son ambos números pequeños. Quizás sea mejor que hagas tres preguntas. Por qué es cierto en teoría, cuando N es prácticamente lo suficientemente grande y si hay algo mal con su código R. Respondí la primera, la segunda es en gran medida calculadora, pero puede haber una buena regla general para este caso específico y la tercera pertenece al desbordamiento de la pila.
Keith
@ Keith lo siento por mis débiles matemáticas, pero ¿puedes mostrar más detalles sobre cómo obtuviste la expectativa?
AdamO