¿Ajustando los datos a las distribuciones?

No soy un estadístico (más de un desarrollador web de investigación) pero he estado escuchando mucho sobre scipy y R en estos días. Así que, por curiosidad, quise hacer esta pregunta (aunque podría parecer una tontería para los expertos de aquí) porque no estoy seguro de los avances en esta área y quiero saber cómo las personas sin antecedentes de estadísticas sólidas abordan estos problemas.

Dado un conjunto de números reales observados en un experimento, digamos que pertenecen a una de las muchas distribuciones que existen (como Weibull, Erlang, Cauchy, Exponential, etc.), ¿existen formas automatizadas de encontrar la distribución correcta y la distribución? parámetros para los datos? ¿Hay algún buen tutorial que me guíe a través del proceso?

Escenario del mundo real: por ejemplo, digamos que inicié una pequeña encuesta y registré información sobre la cantidad de personas con las que una persona habla todos los días, por ejemplo, 300 personas y tengo la siguiente información:

1 10 2 5 3 20 ... ... 

donde XY me dice que esa persona X habló con Y personas durante el período de la encuesta. Ahora, utilizando la información de las 300 personas, quiero encajar esto en un modelo. La pregunta se reduce a: ¿existen formas automatizadas de encontrar los parámetros de distribución y distribución correctos para estos datos o, de no ser así, hay un buen procedimiento paso a paso para lograr lo mismo?

Esta es una pregunta complicada, y no hay respuestas perfectas. Trataré de ofrecerle una visión general de los conceptos principales y señalarle la dirección de algunas lecturas útiles sobre el tema.

Suponga que usted es un conjunto de datos unidimensional y tiene un conjunto finito de funciones de distribución de probabilidad de las que cree que se han generado los datos. Puede considerar cada distribución de forma independiente y tratar de encontrar parámetros que sean razonables dados sus datos. Existen dos métodos para establecer parámetros para una función de distribución de probabilidad a partir de los datos:

  1. Mínimos cuadrados
  2. Máxima verosimilitud

En mi experiencia, la máxima probabilidad se ha preferido en los últimos años, aunque este puede no ser el caso en todos los campos.

Este es un ejemplo concreto de cómo estimar parámetros en R. Considere un conjunto de puntos aleatorios generados a partir de una distribución Gaussiana con una media de 0 y una desviación estándar de 1:

 x = rnorm( n = 100, mean = 0, sd = 1 ) 

Suponga que sabe que los datos se generaron mediante un proceso gaussiano, pero que ha olvidado (¡o nunca supo!) Los parámetros para el gaussiano. Le gustaría usar los datos para obtener estimaciones razonables de la media y la desviación estándar. En R, hay una biblioteca estándar que lo hace muy sencillo:

 library(MASS) params = fitdistr( x, "normal" ) print( params ) 

Esto me dio el siguiente resultado:

  mean sd -0.17922360 1.01636446 ( 0.10163645) ( 0.07186782) 

Esos están bastante cerca de la respuesta correcta, y los números entre paréntesis son intervalos de confianza en torno a los parámetros. Recuerde que cada vez que genere un nuevo conjunto de puntos, obtendrá una nueva respuesta para las estimaciones.

Matemáticamente, esto está usando la máxima probabilidad para estimar tanto la media como la desviación estándar del gaussiano. Probabilidad significa (en este caso) “probabilidad de datos dados los valores de los parámetros”. Máxima probabilidad significa “los valores de los parámetros que maximizan la probabilidad de generar mis datos de entrada”. La estimación de probabilidad máxima es el algoritmo para encontrar los valores de los parámetros que maximizan la probabilidad de generar los datos de entrada, y para algunas distribuciones puede incluir algoritmos de optimización numérica . En R, la mayor parte del trabajo lo realiza fitdistr , que en ciertos casos se denomina óptimo .

Puede extraer la probabilidad de registro de sus parámetros de esta manera:

 print( params$loglik ) [1] -139.5772 

Es más común trabajar con la probabilidad de registro en lugar de evitar errores de redondeo. La estimación de la probabilidad conjunta de sus datos implica multiplicar las probabilidades, que son todas menores que 1. Incluso para un pequeño conjunto de datos, la probabilidad conjunta se aproxima a 0 muy rápidamente, y agregar las probabilidades de registro de sus datos es equivalente a multiplicar las probabilidades. La probabilidad se maximiza a medida que la probabilidad de registro se aproxima a 0 y, por lo tanto, más números negativos se ajustan peor a sus datos.

Con herramientas computacionales como esta, es fácil estimar parámetros para cualquier distribución. Considera este ejemplo:

 x = x[ x >= 0 ] distributions = c("normal","exponential") for ( dist in distributions ) { print( paste( "fitting parameters for ", dist ) ) params = fitdistr( x, dist ) print( params ) print( summary( params ) ) print( params$loglik ) } 

La distribución exponencial no genera números negativos, así que los eliminé en la primera línea. La salida (que es estocástica) se veía así:

 [1] "fitting parameters for normal" mean sd 0.72021836 0.54079027 (0.07647929) (0.05407903) Length Class Mode estimate 2 -none- numeric sd 2 -none- numeric n 1 -none- numeric loglik 1 -none- numeric [1] -40.21074 [1] "fitting parameters for exponential" rate 1.388468 (0.196359) Length Class Mode estimate 1 -none- numeric sd 1 -none- numeric n 1 -none- numeric loglik 1 -none- numeric [1] -33.58996 

La distribución exponencial es en realidad un poco más probable que haya generado estos datos que la distribución normal, probablemente porque la distribución exponencial no tiene que asignar ninguna densidad de probabilidad a números negativos.

Todos estos problemas de estimación empeoran cuando intenta ajustar sus datos a más distribuciones. Las distribuciones con más parámetros son más flexibles, por lo que se ajustarán mejor a sus datos que las distribuciones con menos parámetros. Además, algunas distribuciones son casos especiales de otras distribuciones (por ejemplo, Exponential es un caso especial de Gamma ). Debido a esto, es muy común utilizar el conocimiento previo para restringir los modelos de su elección a un subconjunto de todos los modelos posibles.

Un truco para evitar algunos problemas en la estimación de parámetros es generar una gran cantidad de datos y dejar algunos de los datos para validación cruzada . Para realizar una validación cruzada del ajuste de los parámetros a los datos, deje parte de los datos fuera de su procedimiento de estimación y luego mida la probabilidad de cada modelo en los datos que no se incluyen.

Eche un vistazo a fitdistrplus ( http://cran.r-project.org/web/packages/fitdistrplus/index.html ).

Un par de cosas rápidas a tener en cuenta:

  • Pruebe la función descdist , que proporciona un gráfico de sesgo en comparación con la curtosis de los datos y también muestra algunas distribuciones comunes.
  • fitdist permite ajustar cualquier distribución que pueda definir en términos de densidad y cdf.
  • Luego puede usar gofstat que calcula las estadísticas de KS y AD que miden la distancia del ajuste a partir de los datos.

Es probable que esto sea un poco más general de lo que necesita, pero podría darle algo para continuar.

Una forma de estimar una función de densidad de probabilidad a partir de datos aleatorios es usar una expansión de Edgeworth o Butterworth. Estas aproximaciones utilizan las propiedades de la función de densidad conocidas como cumulantes (los estimadores no sesgados para los cuales son las estadísticas k ) y expresan la función de densidad como una perturbación de una distribución gaussiana.

Ambos tienen algunas debilidades bastante graves como producir funciones de densidad divergentes, o incluso funciones de densidad que son negativas en algunas regiones. Sin embargo, algunas personas los encuentran útiles para datos altamente agrupados, o como puntos de partida para una estimación adicional, o para funciones de densidad estimadas por partes, o como parte de una heurística.

MG Kendall y A. Stuart, La teoría avanzada de la estadística, vol. 1, Charles Griffin, 1963, fue la referencia más completa que encontré para esto, con una enorme página dedicada al tema; la mayoría de los otros textos tenían una oración como máximo o enumeraban la expansión en términos de los momentos en lugar de los cumulantes, lo cual es un poco inútil. Sin embargo, buena suerte al encontrar una copia, tuve que enviar a mi bibliotecario de la universidad a un archivo para que lo archivara … pero esto fue hace años, así que tal vez hoy Internet sea más útil.

La forma más general de su pregunta es el tema de un campo conocido como estimación de densidad no paramétrica , donde se proporciona:

  • datos de un proceso aleatorio con una distribución desconocida, y
  • restricciones en el proceso subyacente

… produce una función de densidad que es más probable que haya producido los datos. (De manera más realista, crea un método para calcular una aproximación a esta función en cualquier punto dado, que puede usar para trabajo adicional, por ejemplo, comparar las funciones de densidad de dos conjuntos de datos aleatorios para ver si podrían provenir de la misma función). proceso).

Personalmente, sin embargo, he tenido poca suerte en el uso de la estimación de la densidad no paramétrica para algo útil, pero si tienes un suministro constante de cordura, deberías estudiarlo.

Básicamente, desea comparar sus datos del mundo real con un conjunto de distribuciones teóricas. Existe la función qqnorm() en la base R, que hará esto para la distribución normal, pero prefiero la función e1071 en e1071 que le permite probar otras distribuciones. Aquí hay un fragmento de código que trazará sus datos reales contra cada una de las distribuciones teóricas que pegamos en la lista. Usamos plyr para revisar la lista, pero también hay otras maneras de hacerlo.

 library("plyr") library("e1071") realData <- rnorm(1000) #Real data is normally distributed distToTest <- list(qnorm = "qnorm", lognormal = "qlnorm", qexp = "qexp") #function to test real data against list of distributions above. Output is a jpeg for each distribution. testDist <- function(x, data){ jpeg(paste(x, ".jpeg", sep = "")) probplot(data, qdist = x) dev.off() } l_ply(distToTest, function(x) testDist(x, realData)) 

No soy un científico, pero si lo hicieras con un lápiz en un papel, la forma obvia sería hacer una gráfica, luego comparar la gráfica con una de una distribución estándar conocida.

Yendo más allá con ese pensamiento, “comparar” es mirar si las curvas de una distribución estándar y las suyas son similares.

Trigonometría, tangentes … sería mi último pensamiento.

No soy un experto, solo otro humilde desarrollador web =)

Para lo que vale, parece que es posible que desee ver la distribución de Poisson.