Estimación de la mediana de la población a partir de un estimador de densidad de núcleo

Question

Tengo una estimación de densidad de núcleo 1-D en forma de dos vectores:
x_grid es un vector de valores de x en los que se muestreó la función de densidad
density es un vector de valores de densidad correspondientes (estimados por una estimación de densidad de núcleo) en cada valor de x en el x_grid.

Si quiero el valor medio de esta función de densidad ($m$ donde $\int_{-\infty}^{m}{f(x) =0.5}$), entonces puedo aproximarlo tomando la integral (usando una técnica numérica como la regla del trapecio) desde $-\infty$ hasta $x_0, x_1, x_2, ...$ hasta que llegue (lo suficientemente cerca) a $0.5$.

Esencialmente, solo estoy sumando las áreas de los trapecios bajo la curva hasta llegar a una suma acumulativa de $0.5$.

¿Cómo puedo hacer algo similar, pero para el valor medio/esperado de la función de densidad?

Answer 1

Aquí hay algunos experimentos utilizando el estimador de densidad de kernel predeterminado en R. Está diseñado de una manera que hace que la estimación de la mediana sea práctica para muestras razonablemente grandes.

He elegido ilustrar con una muestra $V_1, V_2, \dots, V_{500}$ de tamaño $n = 500$ de $\mathsf{Gamma}(\text{forma}=5, \text{tasa}=0.1),$ una distribución sesgada a la derecha. Así que la media poblacional es $\mu = 5/0.1 = 50,$ la moda poblacional es $\delta = (5-1)/0.1 = 40,$ y la mediana poblacional (desde R) es $\eta = 46.71.$ Es factible obtener estimaciones razonables de las tres medidas de centralidad utilizando el KDE predeterminado en R. [Si es necesario, puede encontrar información sobre las distribuciones gamma en Wikipedia.]

qgamma(.5, 5, .1)
[1] 46.70909

Explorando la salida de KDE: La función density(v) produce dos vectores: uno con valores de $x$ uniformemente espaciados y el otro con alturas correspondientes ($y$-valores).

set.seed(418)  # para reproducibilidad
v = rgamma(500, 5, .1)
density(v)

Llamada:
        density.default(x = v)

Datos: v (500 obs.);     Ancho de banda 'bw' = 5.43

       x                 y            
 Min.   : -7.166   Min.   :1.675e-06  
 1st Qu.: 33.912   1st Qu.:2.219e-04  
 Mediana : 74.990   Mediana :2.524e-03  
 Media   : 74.990   Media   :6.080e-03  
 3rd Qu.:116.068   3rd Qu.:1.230e-02  
 Max.   :157.146   Max.   :2.018e-02  

Aquí hay un histograma de densidad de la muestra, junto con el KDE.

lbl = "Densidad de GAMMA(5,0.1) con KDE para n = 500"
hist(v, prob=T, col="skyblue2", main=lbl);  rug(v)
  lines(density(v), lwd=2)

Es conveniente dar nombres a los dos vectores. Además, verificamos que los puntos a lo largo de la cuadrícula horizontal estén igualmente espaciados (llamamos a la distancia común wd), lo que facilita la integración numérica, permitiéndonos ver que el área bajo la curva de KDE es esencialmente $1.$

ht = density(v)$y;  x = density(v)$x
dx = diff(density(v)$x); summary(dx)
   Min. 1st Qu.  Mediana    Media 3rd Qu.    Max. 
 0.3215  0.3215  0.3215  0.3215  0.3215  0.3215 
wd = mean(dx) # solo si todos los 'dx' son iguales
sum(ht*wd)
[1] 1.00097

Medias de muestra, población y de KDE: La media de la muestra $(49.68)$ se calcula a partir de los datos, la media poblacional $(50)$ por fórmula a partir de los parámetros de distribución, y la media de KDE $(49.73)$ mediante integración numérica utilizando el ancho común wd, alturas ht, y valores de $x$ de KDE.

mean(v); 5/.1
[1] 49.68172
[1] 50

 sum(x*ht*wd)
[1] 49.72998

Modos de población y KDE: Por fórmula, el modo poblacional es $4/0.1 = 40.$ El intervalo modal del histograma es $(40,50].$ Y encontramos que el máximo de la curva KDE es aproximadamente $43.96.$ (Con una semilla diferente, la muestra podría haber producido un modo de KDE que no caiga en el intervalo modal del histograma).

4/0.1
[1] 40
mean(x[ht==max(ht)])
[1] 43.96074

Medianas de muestra, población y de KDE: La mediana de la muestra $(46.94)$ se encuentra directamente a partir de la muestra con median, la mediana poblacional $(46.71)$ mediante integración numérica con qgamma, y la mediana de KDE $(47.50)$ mediante integración numérica (hasta la mediana) a partir de la salida de KDE. Los tres valores están cerca de $47.$

median(v);  qgamma(.5, 5, .1)
[1] 46.93721
[1] 46.70909
kde.med = min(x[cumsum(ht*wd)>=.5]);  kde.med
[1] 47.49778

Funciones de distribución empírica, poblacional y de KDE: La ECDF de una muestra de tamaño $n$ tiene saltos de tamaño $1/n$ en los valores ordenados de la muestra (con saltos de $k/n$ en datos discretos o redondeados $k$ observaciones empatadas en el mismo valor).

La figura superior a continuación compara la ECDF con la curva de CDF poblacional; las líneas sólidas muestran la ubicación de la mediana poblacional, y las líneas punteadas la ubicación de la mediana de la muestra. De manera similar, la figura inferior compara la CDF de KDE y la CDF poblacional (y las respectivas medianas). En caso de interés, se muestra a continuación mi código R para la figura.

par(mfrow=c(2,1))  # permitir dos gráficos por panel
 hdr1 = "ECDF de Muestra (Puntos) y CDF de GAMMA(5,.1)"
 plot(ecdf(v), main=hdr1)
    lines(c(-10,median(v), median(v)), c(.5, .5, 0), lwd=2, lty="dashed")
  curve(pgamma(x, 5, .1), add=T, lwd=3, col="green2")
    lines(c(-10, qgamma(.5, 5, .1), qgamma(.5, 5, .1)), c(.5, .5, 0), col="red")
 hdr2 = "CDF de KDE y CDF de GAMMA(.5, .1)"
 plot(x, kde.cdf, pch=20, main=hdr2)
     lines(c(-10, kde.med, kde.med), c(.5,.5,0), lwd=2, lty="dashed")
   curve(pgamma(x, 5, .1), add=T, lwd=2, col="green2")
     lines(c(-10, qgamma(.5, 5, .1), qgamma(.5, 5, .1)), c(.5, .5, 0), col="red")
par(mfrow=c(1,1))  # volver a la representación de figura única

Notas: (1) Para su KDE, necesita determinar si los puntos de la cuadrícula horizontal están igualmente espaciados como en este ejemplo con el KDE predeterminado en R. Si no tendrá que usar anchos de rectángulo variables en la integración numérica. (2) No restringí el KDE a valores positivos. Es una buena idea restringir el KDE para que coincida con el soporte de la distribución involucrada. (3) Tenga en cuenta que los resúmenes de los valores de $x$ y $y$ dados en la primera sección de esta Respuesta no son directamente aplicables para estimar la media, mediana y moda de la población. (4) Descargo de responsabilidad: Gran parte de lo que sé sobre los detalles particulares del KDE predeterminado de R fue aprendido haciendo esta investigación. Las personas en este sitio que saben más sobre el diseño de la función density en R, pueden tener comentarios importantes para añadir.