Estimador robusto de los parámetros de una distribución de mezcla

Question

Tengo una variable aleatoria que toma uno de tres valores en una progresión aritmética, $a$ , $a+b$ , $a+2b$ , con probabilidades $1/4,1/2,1/4$ , donde $b\approx 0.5a$ a $2a$ . Está perturbado por un poco de ruido (digamos gaussiano con desviación igual a una décima de $a$ ), y puede haber una pequeña proporción de valores atípicos, por ejemplo $10\%$ .

Me gustaría construir un estimador robusto de $a$ et $b$ de una muestra de $N$ valores, donde $N$ es del orden de $50$ .

Podría, por ejemplo, intentar encontrar tres modos en el histograma, y luego clasificar los valores asignándolos al modo más cercano y estimar las medias de las clases. También supongo que la mediana sería un buen estimador de $a+b$ . (Y la distribución de $|v-\overline{a+b}|$ sería bimodal).

¿Puede aconsejar un enfoque para $b$ ?

---

Siguiendo los comentarios de @whuber, puedo reformular la pregunta como la estimación de las medias de una mezcla gaussiana de tres miembros, sabiendo que estas medias son equidistantes, las varianzas son iguales y hay una buena separación entre los miembros (pero algunos valores atípicos).

Answer 1

Un análisis de conglomerados, seguido de la búsqueda de las medianas de los conglomerados, probablemente funcionaría bien. Sin embargo, existen procedimientos sencillos, eficaces y demostrables que podrían recomendarse porque (1) son relativamente fáciles de estudiar y analizar; (2) aprovechan la progresión aritmética asumida y la suposición de que aproximadamente la mitad de los datos se encuentran cerca de la modalidad media; y (3) pueden ajustarse para lograr la robustez deseada (dentro de unos límites).

Un procedimiento comienza por estimar $a+b$ como la mediana de los datos. Entonces se espera que los residuos absolutos estén agrupados, en proporciones aproximadamente iguales, cerca de $0$ et $b$ . Especificando un cuantil adecuado $q$ entre $0.5$ et $1.0$ puede obtener un presupuesto preliminar $b^{*}$ de $b$ Por ejemplo, el $0.75$ cuantil debería ser una buena estimación. Con esta estimación preliminar en la mano, complete un simple análisis de conglomerados dividiendo los residuos absolutos en aquellos menores de $b^{*}$ y los mayores de $b^{*}$ . La mediana de este último grupo es una estimación robusta de $b.$ La estimación de $a$ se obtiene restando la estimación de $b$ a partir de la estimación de $a+b.$ (Ver el código al final de este post).

Hay una compensación: si $q$ no es lo suficientemente alto, puede haber una posibilidad apreciable de que este cuantil se sitúe en el grupo cercano a $0$ y, por lo tanto, subestimar groseramente $b$ . Si $q$ es demasiado alto, el cuantil será sensible a los valores atípicos y quizás sobreestimará groseramente $b.$ Ajustar $q$ para lograr el equilibrio deseado entre estos comportamientos permite afinar el procedimiento.

He aquí un análisis rápido. La posibilidad de que $b^{*}$ es demasiado pequeño en un conjunto de datos de tamaño $n$ es aproximadamente la probabilidad de que más de $qn$ de los datos provienen del componente situado en $a+b$ . Ignorando los valores atípicos -la mayoría de los cuales probablemente no se acerquen a $a+b,$ de todos modos esto tiene un Binomio $(n, 1/2)$ distribución. Por ejemplo, con $n=50$ et $q=0.75$ la posibilidad es sólo $0.00015.$ Eso puede ser lo suficientemente apreciable (me pasó una vez en un conjunto de $200$ conjuntos de datos simulados) que podría querer aumentar $q$ a algún lugar entre $0.8$ et $0.9.$ Sin embargo, esto reducirá el punto de ruptura: con $q=0.9,$ por ejemplo, más de cinco valores atípicos pueden cambiar radicalmente $b^{*}.$

He realizado una simulación con $n=50$ et $q=0.875,$ utilizando perturbaciones iguales a $0.1b$ veces una distribución de Cauchy. Es una prueba bastante severa porque esta distribución produce valores atípicos a un gran ritmo: tiene alrededor de un $1$ en $8$ posibilidad de cambiar un valor en más de $b/2,$ acercándolo así a otro modo (o alejándolo de todos los modos).

Aquí están los histogramas de la primera $20$ iteraciones. He arreglado $a=b=1$ y han marcado los lugares $a=1,a+b=2,a+2b=3$ con guiones negros verticales. Las estimaciones se muestran con líneas rojas. Están en buena concordancia con los valores reales para estas iteraciones. (Las elecciones de $a$ et $b$ no importan: lo único relevante es la dispersión del ruido aleatorio, expresada en una escala donde $b$ es una unidad).

Los valores atípicos son evidentes en las iteraciones 5, 12 y 14 porque son muy extremos. También se pueden ver muchos datos que caen vagamente entre dos de las modalidades, lo que dificulta su clasificación.

La siguiente figura es un gráfico de dispersión para mostrar la relación entre los valores reales de los modos (es decir, $1,2,3$ ) y los valores estimados de las 2000 iteraciones.

Evidentemente ninguna de las estimaciones era mala y suelen ser imparciales. La gran mayoría estaban bien dentro de $0.1$ de los valores reales. (Sobre $92\%$ de las estimaciones de $a$ estaban entre $0.9$ et $1.1$ mientras que casi $99\%$ de las estimaciones de $b$ estaban entre $0.9$ et $1.1$ .)

Por último, veamos cómo las estimaciones de $a$ et $b$ están relacionados, dibujando su diagrama de dispersión:

Como es de esperar, están correlacionados negativamente. (La curva es una suavidad local no lineal, por lo que llama la atención que sea casi lineal en todo el rango de estimaciones). Como la nube se concentra en $(1,1)=(a,b),$ las estimaciones son insesgadas. (Sus medias en esta simulación fueron $1.009$ et $0.993,$ respectivamente, con un coeficiente de correlación de $-60\%.$ )

Finalmente, el R El código para crear estas estimaciones es rápido y sencillo. Devuelve las estimaciones de los tres modos. Se garantiza que la estimación media está exactamente a medio camino entre las dos estimaciones extremas. (Esto difiere de un procedimiento general para estimar las ubicaciones de una mezcla de tres componentes. También difiere en que se basa en la suposición de que el peso del componente medio es la mitad del total).

estimator <- function(x, q=0.875) {      # `x` is an array of data
  ab.hat <- median(x)                    # Estimate of a+b
  z <- abs(x - ab.hat)                   # Absolute residuals
  b.hat <- quantile(z, q)                # Preliminary estimate of b
  threshold <- b.hat / 2                 # Cutoff between 0 and b
  b.hat <- median(z[z >= threshold])     # Final estimate of b
  c(ab.hat-b.hat, ab.hat, ab.hat+b.hat)
}

Answer 2

Propuesta rápida. Siéntase libre de mejorarla.

Empiece por su solución parcial, es decir, deje que $m_2$ sea su mediana general. Como sabe que aproximadamente la mitad de los puntos deben estar en el grupo medio, eliminar la mitad de los puntos, los más cercanos a su mediana general, podría ser una buena idea.

Llegados a este punto, te quedarás con la otra mitad, que deberás dividir en dos grupos de aproximadamente el mismo tamaño. Tomando la mediana de todos los puntos que son menores que $m_2$ , nómbralo $m_1$ y luego la mediana de todos los puntos que son mayores que $m_2$ , nómbralo $m_3$ es una opción interesante.

Por último, dejemos que $(a,b) = \mathrm{argmin} \sum_{i=1}^3 \{a+b(i-1) - m_i\}^2$ .

Esta solución puede no ser la mejor si sus grupos no poseen variantes similares, o variantes enormes de manera que se solapen mucho.