[R-es] Goodness

Carlos J. Gil Bellosta cgb en datanalytics.com
Mie Feb 10 22:59:10 CET 2010


Hola, ¿qué tal?

En efecto, es un tema de discusión recurrente y creo que fui yo quien lo
saqué a colación por primera vez el año pasado. Trabajaba en un lugar en
el que tenían un programa (en Matlab, aunque lo que se aplique a Matlab
se puede extrapolar a R) que trataba de resolver un problema análogo:
dada una muestra de valores (eran siempre positivos) buscar la
distribución (dentro de un conjunto preestablecido) que mejor se ajusta
a ellos.

Además, en realidad, no interesaba conocer la distribución sino para
calcular posteriormente ciertos estadísticos asociados a ella. Se
trataba de determinados cuantiles.

La clase de distribuciones que manejábamos incluía algunas para las que
existían parámetros que se podían calcular analíticamente. Otros no y
hacía falta recurrir a el equivalente en Matlab de "fitdist" (que llama
internamente a "optim", si mal no recuerdo).

Además, usar fitdist "universalmente" daba lugar a respuestas
manifiestamente subóptimas (caso eminente: distribuciones para las que
uno de los parámetros indicaba el soporte de la distribución). Mis
predecesores en el desarrollo de esa solución tuvieron el acierto (que
les evitó mucho trabajo) de suprimir los "warnings" para no tener que
preocuparse de los alarmantes avisos de Matlab y así llegar a casa a
cenar a tiempo todos los días.

Después de algunas discusiones (entre ellas una en la lista), la
política que se siguió para el desarrollo de la solución fue el siguiente:

1) Para el ajuste de las distribuciones, crear clases que las extendían
y que incorporaban una nuevo método (se llamaba "mle") que era capaz de
identificar aquellos parámetros calculables analíticamente y que
realizaba llamadas a "optim" para los restantes. Cada distribución, por
lo tanto, se ajustaba, potencialmente, de una manera distinta.

2) Para solventar el problema del posible sobreajuste (y eludir la
falacia del polinomio de orden n-1) y dado que el interés máximo
radicaba en determinados cuantiles, lo que se ensayó fue realizar varios
ajustes con muestras de los datos originales, calcular el estadístico de
interés para cada una de ellas y utilizar la varianza de la muestra (de
estadísticos) resultante como índice de sobreajuste. La combinación de
los estadísticos tradicionales de ajuste con el de la varianza del
estadístico resultante guiaba la elección de la distribución "ganadora".

El problema de encontrar la "mejor" distribución, aunque escandalice a
los puristas, aparece tantas veces en la práctica que no cabe sino
afrontarlo. Y cada vez, en cada contexto, con una estrategia distinta.

Un cordial saludo,

Carlos J. Gil Bellosta
http://www.datanalytics.com





Pablo Emilio Verde wrote:
> Hola,
> 
> Este es un tema con discusion recurrente. El año pasado mande dos scripts
> en R, uno que seleccion el model con minimo AIC dentro de una familia de
> distribuciones.
> 
> El otro script usa simulacion para validar modelos y omite el uso de test
> de bondad de ajuste (que son obsoletos).
> 
> Ruben, fijate en los archivos del año pasado y seguro encontras estos
> scripts, quizas te sirvan.
> 
> Saludos,
> 
> Pablo
> 
> 
> ----- Original Message -----
> From: "J. Miguel Marin" <jmmarin en est-econ.uc3m.es>
> To: <r-help-es en r-project.org>
> Sent: Tuesday, February 09, 2010 7:19 PM
> Subject: Re: [R-es] Goodness
> 
> 
> Hola,
> leyendo entre líneas yo creo que Daniel busca algo parecido a lo que
> tiene Statgraphics, que tiene una batería de distribuciones para
> ajustar los datos "a tu gusto".
> Lo más parecido que tiene R es justamente el "fitdistr".
> En realidad, yo creo que Daniel busca una generalización de un test de
> Kolmogorov-Smirnov "para cualquier cosa", aunque me temo que a pesar de
> que eso es efectivamente erróneo mucha gente de Ciencias lo usa en plan
> ciego.
> 
> 
> 
>> Estoy de acuerdo con Rubén.
>>
>> El planteamiento de Daniel me recuerda a una pregunta típica en
>> Médicos, Biólogos, Ingenieros, con poca experiencia en Estadística.
>> La pregunta a menudo es: ¿Qué modelo de regresión no lineal será el
>> que mejor se ajusta para explicar mis observaciones de Y en función
>> de X? y con algo de sorna siempre contesto: "El polinomio con n-1
>> datos".
>>
>> Creo que este caso es equivalente: La distribución que mejor se
>> ajuste a una colección de datos no le veo sentido. Suponiendo que se
>> han obtenido mediante muestreo, se puede dar un ajuste espúreo a
>> cualquier distribución. Algunas convergen bastante bien sobre otras
>> bajo determinadas condiciones.
>>
>> En una distribución creo que es preciso plantear el modelo de
>> generación de los datos o si ello no es posible conformarse con una
>> estimación de la densidad mediante procedimientos como las
>> estimaciones basadas en funciones núcleo. En concreto la función
>> "density" del paquete stats estima la densidad de una variable a
>> partir de un conjunto de datos mediante algunas de las funciones
>> núcleo más habituales.
>>
>> Proponer una función de distribución con toda su parsimonia en los
>> parámetros sin justificación en la generación de los datos a mí me
>> parece poco defendible a partir de una mera falta de desajuste de las
>> observaciones respecto de un modelo del que no se justifica por qué
>> este y no otro.
>>
>> Si las observaciones proceden de una mezcla de distribuciones, no
>> habrá ningún modelo sencillo, del grupo de los que podríamos llamar
>> básicos, que realmente prediga las observaciones y no le veo ningún
>> sentido a la propuesta de Daniel.
>>
>>
>> Rubén Roa escribió:
>>>  -----Mensaje original-----
>>> De: r-help-es-bounces en r-project.org
>>> [mailto:r-help-es-bounces en r-project.org] En nombre de daniel pacheco
>>> gomez Enviado el: martes, 09 de febrero de 2010 12:02
>>> Para: r-help-es en r-project.org
>>> Asunto: [R-es] Goodness
>>>
>>>
>>> Hola,
>>>
>>> LLevo buscando desde hace tiempo como hacer el Goodness of fit test
>>> en R. Es decir, me explico, intento hacer una cosa parecida que se
>>> hace en Minitab, por ejemplo, yo tengo un conjunto de datos, y  lo
>>> que quiero es sabes que tipo de distibución es, en minitab se hace
>>> un histograma para ver si se ajusta bien o no a la campana de Gauss,
>>> luego vemos si aproximar la distribución de la muestra por una
>>> Normal es lícito. Es decir, minitab genera un Probability Plot of
>>> Adjusted Value y devuelve el Anderson Darling y el P-Value, al hacer
>>> esto se observa que el p-value del contraste de hipótesis utilizado
>>> es menor que el nivel de significación estandar(0,05), luego no
>>> podemos aceptar que la hipótesis  nula de que la muestra provenga de
>>> una distribución normal.
>>>
>>> Es entonces cuando en minitab se le aplica el "Goodness of Fit
>>> Test", y se observe que distribución se ajusta mejor a la muestra
>>> obtenida. Entonces se le aplica el Process capability según la
>>> distribución y se analizan los resultados.
>>>
>>> Mi duda es, ¿exite alguna función-libreria-paquete, que al pasarle
>>> la muestra, te diga que tipo de distribución es?.
>>>
>>> ---
>>> No lo creo, no es está en la filosofía de R una forma tan automática
>>> de responder tus preguntas.
>>> Lo que tú quieres hacer yo lo hago con fitdistr de MASS, y si tengo
>>> dos o mas distribuciones que pueden explicar mis datos, uso el AIC.
>>> Lo de qué distribución corresponde (nota la diferencia con 'mejor
>>> ajusta') a tus datos, lo puedes elucidar pensando en la naturaleza
>>> de tus datos.
>>> Son continuos o discretos?, univariados o multivariados?.
>>> Si son discretos, son conteos? (Poisson?, binomial negativa?), son
>>> binarios? (binomial)?, tienen exceso de ceros? (sobredispersos)?
>>> Si son continuos, admiten cero y valores negativos? (normal?), sólo
>>> valores positivos? (lognormal {proceso multiplicativo}, gamma {suma
>>> de exponenciales}?), son distancias o tiempos entre eventos?
>>> (exponencial?), etc.
>>>
>>>
>>> par(mfrow=c(2,1),oma=c(2,2,1,1),mar=c(2,2,1,1))
>>>
>>> #Ejemplo 1 (solución analítica para los emv):
>>>
>>> x <- rnorm(250, 5, 3)
>>> hist(x, prob=TRUE, ylim=c(0,1.25*max(hist(x,plot=FALSE)$density)))
>>> curve(dnorm(x,4.5,3.2,),col="blue",add=TRUE)
>>> library(MASS)
>>> x.likfit <- fitdistr(x,'normal') #sin proporcionar valores iniciales
>>>
> curve(dnorm(x,mean=x.likfit$estimate[1],sd=x.likfit$estimate[2],),col="red",
> add=TRUE)
>>> #Ejemplo 2 (solución numérica para los emv):
>>>
>>> x <- rgamma(250, 2, 25)
>>> hist(x, prob=TRUE, ylim=c(0,1.25*max(hist(x,plot=FALSE)$density)))
>>> curve(dgamma(x,shape=3,rate=30),col="blue",add=TRUE)
>>> library(MASS)
>>> x.likfit <- fitdistr(x,'gamma',start=list(shape =3, rate=30)) #con
>>> valores iniciales
>>>
> curve(dgamma(x,shape=x.likfit$estimate[1],rate=x.likfit$estimate[2],),col="r
> ed",add=TRUE)
>>> HTH
>>> Rubén
>>>
>>>
> ____________________________________________________________________________
> ________ Dr. Rubén
>>> Roa-Ureta
>>> AZTI - Tecnalia / Marine Research Unit
>>> Txatxarramendi Ugartea z/g
>>> 48395 Sukarrieta (Bizkaia)
>>> SPAIN
>>>
>>> _______________________________________________
>>> R-help-es mailing list
>>> R-help-es en r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>>
>> _______________________________________________
>> R-help-es mailing list
>> R-help-es en r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>
> 
> 
> 
> 
> jm~
> 
> _______________________________
> 
>         J. Miguel Marin
> 
> http://www.est.uc3m.es/jmmarin
> 
>     Dep. of Statistics
> University Carlos III of Madrid
>         Spain (E.U.)
> 
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
> 
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>



Más información sobre la lista de distribución R-help-es