¿Cómo encontrar el modo de una función de densidad de probabilidad?

Inspirado por mi otra pregunta , me gustaría preguntar cómo se encuentra el modo de una función de densidad de probabilidad (PDF) de una función $f(x)$ .

¿Hay algún procedimiento de "libro de cocina" para esto? Aparentemente, esta tarea es mucho más difícil de lo que parece al principio.

Decir "el modo" implica que la distribución tiene uno y solo uno. En general, una distribución puede tener muchos modos, o (posiblemente) ninguno.

Si hay más de un modo, debe especificar si desea todos ellos o solo el modo global (si hay exactamente uno).

Suponiendo que nos limitemos a distribuciones unimodales *, para que podamos hablar del modo "the", se encuentran de la misma manera que encontrar máximas de funciones de manera más general.

* tenga en cuenta que la página dice " como el término" modo "tiene múltiples significados, también lo hace el término" unimodal " " y ofrece varias definiciones de modo, que pueden cambiar lo que, exactamente, cuenta como un modo, ya sea que haya 0 1 o más, y también altera la estrategia para identificarlos. Observe particularmente cuán general es la redacción "más general" de lo que es la unimodalidad en el párrafo inicial " unimodalidad significa que solo hay un valor más alto, definido de alguna manera "

Una definición ofrecida en esa página es:

Un modo de distribución de probabilidad continua es un valor en el que la función de densidad de probabilidad (pdf) alcanza su valor máximo

Así dado una definición específica del modo, la encontrará como encontraría esa definición particular de "valor más alto" cuando se trata de funciones de manera más general (suponiendo que la distribución es unimodal bajo esa definición).

Hay una variedad de estrategias en matemáticas para identificar tales cosas, dependiendo de las circunstancias. Vea, la sección "Encontrar máximos y mínimos funcionales" de la página de Wikipedia sobre Máximos y mínimos que ofrece una breve discusión.

Por ejemplo, si las cosas son lo suficientemente buenas, digamos que estamos tratando con una variable aleatoria continua, donde la función de densidad tiene una primera derivada continua, puede proceder tratando de encontrar dónde la derivada de la función de densidad es cero, y verificando de qué tipo de punto crítico es (máximo, mínimo, punto de inflexión horizontal). Si hay exactamente uno de esos puntos que es un máximo local, debería ser el modo de una distribución unimodal.

Sin embargo, en general, las cosas son más complicadas (por ejemplo, el modo puede no ser un punto crítico), y entran en juego las estrategias más amplias para encontrar el máximo de funciones.

A veces, encontrar dónde los derivados son cero algebraicamente puede ser difícil o al menos engorroso, pero aún puede ser posible identificar máximos de otras maneras. Por ejemplo, puede ser que uno invoque consideraciones de simetría al identificar el modo de una distribución unimodal. O uno podría invocar alguna forma de algoritmo numérico en una computadora, para encontrar un modo numérico.

Aquí hay algunos casos que ilustran cosas típicas que debe verificar, incluso cuando la función es unimodal y al menos continua por partes.

Entonces, por ejemplo, debemos verificar los puntos finales (diagrama central), los puntos donde la derivada cambia de signo (pero puede no ser cero; primer diagrama) y los puntos de discontinuidad (tercer diagrama).

En algunos casos, las cosas pueden no ser tan ordenadas como estas tres; debe tratar de comprender las características de la función particular con la que está tratando.

No he tocado el caso multivariante, donde incluso cuando las funciones son bastante "agradables", solo encontrar máximos locales puede ser sustancialmente más complejo (por ejemplo, los métodos numéricos para hacerlo pueden fallar en un sentido práctico, incluso cuando lógicamente deben tener éxito finalmente).

Esta respuesta se centra completamente en la estimación de modo de una muestra, con énfasis en un método en particular. Si hay algún sentido fuerte en el que ya conoces la densidad, analítica o numéricamente, entonces la respuesta preferida es, en resumen, buscar el máximo máximo único o múltiple directamente, como en la respuesta de @Glen_b.

Los "modos de media muestra" pueden calcularse utilizando la selección recursiva de la media muestra con la longitud más corta. Aunque tiene raíces más largas, Bickel y Frühwirth (2006) presentaron una excelente presentación de esta idea.

La idea de estimar el modo como el punto medio del intervalo más corto que contiene un número fijo de observaciones se remonta al menos a Dalenius (1965). Ver también Robertson y Cryer (1974), Bickel (2002) y Bickel y Frühwirth (2006) sobre otros estimadores del modo.

Las estadísticas de orden de una muestra de valores de x están definidas por x ( 1 ) ≤ x ( 2 ) ≤ ⋯ ≤ x ( n - 1 ) ≤ x ( n $n$$x$$x_{(1)} \le x_{(2)} \le \cdots \le x_{(n-1)} \le x_{(n)}$ .

El modo de media muestra se define aquí usando dos reglas.

Regla 1. Si , el modo de media muestra es x ( 1 ) . Si n = 2 , el modo de media muestra es ( x ( 1 ) + x ( 2 ) ) / 2 . Si n = 3 , el modo de medio-muestra es ( x ( 1 ) + x ( 2 ) ) / 2 si x ( 1 ) y x ( $n = 1$$x_{(1)}$$n = 2$$(x_{(1)} + x_{(2)}) / 2$$n = 3$$(x_{(1)} + x_{(2)}) / 2$$x_{(1)}$$x_{(2)}$están más cerca que y x ( 3 ) , ( x ( 2 ) + x ( 3 ) ) / 2 si lo opuesto es verdad, y x ( 2 $x_{(2)}$$x_{(3)}$$(x_{(2)} + x_{(3)}) / 2$$x_{(2)}$ de otro modo.

Regla 2. Si , aplicamos la selección recursiva hasta que quede con 3 o menos valores. Primero deje h 1 = ⌊ n / 2 ⌋ . La mitad más corta de los datos del rango k al rango k + h 1 se identifica para minimizar x ( k + h 1 ) - x ( k ) sobre k = 1 , ⋯ , n - h $n \ge 4$$3$$h_1 = \lfloor n / 2\rfloor$$k$$k + h_1$$x_{(k + h_1)} - x_{(k)}$$k = 1, \cdots, n - h_1$ . Entonces la mitad más corta de esas valores se identifica usando h 2 = ⌊ h 1 / 2 ⌋ , y así sucesivamente. Para terminar, use la Regla 1.$h_1 + 1$$h_2 = \lfloor h_1 / 2\rfloor$

La idea de identificar la mitad más corta se aplica en el "shorth" nombrado por JW Tukey e introducido en el estudio de robustez Princeton de estimadores de ubicación por Andrews, Bickel, Hampel, Huber, Rogers y Tukey (1972, p.26) como el media de la media longitud más corta para h = ⌊ n / 2 ⌋ . Rousseeuw (1984), basándose en una sugerencia de Hampel (1975), señaló que el punto medio de la mitad más corta ( x k + x ( k + h $x_{(k)}, \cdots, x_{(k + h)}$$h = \lfloor n / 2 \rfloor$ es el estimador de ubicación de mediana menor de cuadrados (LMS) para x . Ver Rousseeuw (1984) y Rousseeuw y Leroy (1987) para aplicaciones de LMS e ideas relacionadas a la regresión y otros problemas. Tenga en cuenta que este punto medio LMS también se llama shorth en algunas publicaciones más recientes (por ejemplo, Maronna, Martin y Yohai 2006, p.48). Además, la mitad más corta en sí misma a veces también se llama shorth, como lo indica el título de Grübel (1988). Para una implementación de Stata y más detalles, ver desde SSC.$(x_k + x_{(k + h)}) / 2$$x$shorth

Algunos comentarios generales siguen las ventajas y desventajas de los modos de media muestra, desde el punto de vista de los analistas de datos prácticos tanto como los estadísticos matemáticos o teóricos. Cualquiera sea el proyecto, siempre será prudente comparar los resultados con medidas de resumen estándar (p. Ej., Medianas o medias, incluidas las medias geométricas y armónicas) y relacionar los resultados con gráficos de distribuciones. Además, si su interés está en la existencia o el alcance de la bimodalidad o la multimodalidad, será mejor mirar directamente a estimaciones adecuadamente suavizadas de la función de densidad.

Estimación del modo Al resumir dónde están los datos más densos, el modo de media muestra agrega un estimador automático del modo a la caja de herramientas. Las estimaciones más tradicionales del modo basadas en la identificación de los picos en los histogramas o incluso los gráficos de densidad del núcleo son sensibles a las decisiones sobre el origen o el ancho del contenedor o el tipo de núcleo y el ancho medio del núcleo y, en cualquier caso, son más difíciles de automatizar. Cuando se aplica a distribuciones que son unimodales y aproximadamente simétricas, el modo de media muestra estará cerca de la media y la mediana, pero más resistente que la media a los valores atípicos en cualquiera de las colas. Cuando se aplica a distribuciones que son unimodales y asimétricas, el modo de media muestra estará más cerca del modo identificado por otros métodos que la media o la mediana.

Simplicidad La idea del modo de media muestra es bastante simple y fácil de explicar a estudiantes e investigadores que no se consideran especialistas en estadística.

Interpretación gráfica El modo de media muestra puede relacionarse fácilmente con las pantallas estándar de distribuciones, como las gráficas de densidad del grano, la distribución acumulativa y las gráficas de cuantiles, histogramas y gráficas de tallo y hojas.

Al mismo tiempo, tenga en cuenta que

No es útil para todas las distribuciones Cuando se aplica a distribuciones que tienen aproximadamente forma de J, el modo de media muestra se aproximará al mínimo de los datos. Cuando se aplica a distribuciones que tienen aproximadamente forma de U, el modo de media muestra estará dentro de la mitad de la distribución que tenga una densidad promedio más alta. Ninguno de los comportamientos parece especialmente interesante o útil, pero igualmente hay poca necesidad de resúmenes de modo único para distribuciones en forma de J o en forma de U. Para las formas en U, la bimodalidad hace que la idea de un solo modo sea discutible, si no es inválida.

Lazos La mitad más corta puede no estar definida de manera única. Incluso con datos medidos, el redondeo de los valores informados con frecuencia puede dar lugar a vínculos. Qué hacer con dos o más mitades más cortas ha sido poco discutido en la literatura. Tenga en cuenta que las mitades atadas pueden superponerse o ser disjuntas.

hsmode$t$$t$$\lceil t/ 2\rceil$

$-9, -4, -1 , 0, -1, 4, 9$$-0.5$$0$$1 + \lfloor n / 2\rfloor$$n$$n$, lo cual es difícil de lograr dados otros objetivos, en particular que la longitud de la ventana nunca debe disminuir con el tamaño de la muestra. Preferimos creer que este es un problema menor con conjuntos de datos de tamaño razonable.

$1 + \lfloor n / 2\rfloor$$n$$n$$n = 1,$$n = 2$$\lceil n / 2\rceil$

$1.6, 3.11, 3.95, 4.2, 4.2, 4.62, 4.62, 4.62, 4.7, 4.87, 5.04, 5.29, 5.3, 5.38, 5.38, 5.38, 5.54, 5.54, 5.63, 5.71, 6.13, 6.38, 6.38, 6.67, 6.69, 6.97, 7.22, 7.72, 7.98, 7.98, 8.74, 8.99, 9.27, 9.74, 10.66.$ The Stata implementation hsmode reports a mode of 5.38. Robertson and Cryer's own estimates using a rather different procedure are $5.00, 5.02, 5.04$. Compare with your favourite density estimation procedure.

Andrews, D.F., P.J. Bickel, F.R. Hampel, P.J. Huber, W.H. Rogers and J.W. Tukey. 1972. Robust estimates of location: survey and advances. Princeton, NJ: Princeton University Press.

Bickel, D.R. 2002. Robust estimators of the mode and skewness of continuous data. Computational Statistics & Data Analysis 39: 153-163.

Bickel, D.R. and R. Frühwirth. 2006. On a fast, robust estimator of the mode: comparisons to other estimators with applications. Computational Statistics & Data Analysis 50: 3500-3530.

Dalenius, T. 1965. The mode - A neglected statistical parameter. Journal, Royal Statistical Society A 128: 110-117.

Grübel, R. 1988. The length of the shorth. Annals of Statistics 16: 619-628.

Hampel, F.R. 1975. Beyond location parameters: robust concepts and methods. Bulletin, International Statistical Institute 46: 375-382.

Maronna, R.A., R.D. Martin and V.J. Yohai. 2006. Robust statistics: theory and methods. Chichester: John Wiley.

Robertson, T. and J.D. Cryer. 1974. An iterative procedure for estimating the mode. Journal, American Statistical Association 69: 1012-1016.

Rousseeuw, P.J. 1984. Least median of squares regression. Journal, American Statistical Association 79: 871-880.

Rousseeuw, P.J. and A.M. Leroy. 1987. Robust regression and outlier detection. New York: John Wiley.

This account is based on documentation for

Cox, N.J. 2007. HSMODE: Stata module to calculate half-sample modes, http://EconPapers.repec.org/RePEc:boc:bocode:s456818.

See also David R. Bickel's website here for information on implementations in other software.

If you have samples from the distribution in a vector "x", I would do:

 mymode <- function(x){
   d<-density(x)
   return(d$x[which(d$y==max(d$y)[1])])
 }

You should tune the density function so it is smooth enough on the top ;-).

If you have only the density of the distribution, I would use an optimiser to find the mode (REML, LBFGS, simplex, etc.)...

 fx <- function(x) {some density equation}
 mode <- optim(inits,fx)

Or use a Monte-Carlo sampler to get some samples from the distribution (package rstan) and use the procedure above. (Anyway, Stan package as an "optimizing" function to get the mode of a distribution).