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?. He probado hacer varias cosas, aunque no me termina de convencer la manera de hacerla, ni los datos obtenidos. 1 - Una es con la libreria "nortest", que la encontre en un documento que se llama "FITTING DISTRIBUTIONS WITH R" Ricci, en ella viene una función que se llama "ad.test()", donde al pasarle la muestra te devuelve un AD(anderson darling) y p-value de la distribución, entonces comparo el p-value con la normal, y si cuadra, pues teoricamente deberia de probar con las otras distribuciones. Es decir, deberia de calcular los parametros para la "rweibull" por ejemplo, y aplicarla por ejemplo a 100 iteraciones, y volver a pasarle el ad.test(), y ver si el p-value obtenido coincide. Aunque esta forma de hacerlo no me convence, pues tendria que ir probando una a una cada distribucion, ir calculandole los parametros y tal, seria una forma muy manual. Tampoco se si para calcular la distribucion que es, deberia aplibarle por ejemplo en la weibul, la rweibul, pweibul.. 2 - Otra forma de hacerlo, o mas o menos un camino seguir, fue con la libreria "qAnalyst", la cual tiene un par de funciones para esto, por ejemplo, la funcion "funInfoFun", donde se le va pasando la muestra de datos, y la distribucion, y creo, que directamente te da los valores de AD y p-value, lo malo es que esta funcion solo es para distribuciones normal, lognormal, gamma, weibull, logistic, lo cual si no son ninguna de estas no sabria que hacer, y que al comparar el resultado con el primer metodo no me coinciden los datos. Un ejemplo seria: data(warpTiles) x=warpTiles$warping infoX=funInfoFun(x,"weibull") andersonDarlingFun(x=x, infoX$densfun, infoX$theta) Lo bueno de esta libreria es para calcular el process capability segun la distribucion, es decir, las funciones "capabilityNormal" y "capabilityNotNormal", la cual esta ultima habria que pasarle el tipo de distribución. Bueno, como veis, me lio bastante con esto, perdonad por el tocho escrito, pero era para ser lo mas detallado y descriptivo posible. Si alguien ha hecho esto alguna vez, sabe alguna otra manera de hacerlo, sabe como hacerlo con alguno de estos dos metodos pero bien, o cualquier cosa, pues le estaria muy agradecio. Asi que gracias de antemano a tod@s. Navega con el navegador más seguro de todos. ¡Descárgatelo ya! _________________________________________________________________ [[elided Hotmail spam]] inero! [[alternative HTML version deleted]]
-----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="red",add=TRUE)
HTH
Rubén
____________________________________________________________________________________
Dr. Rubén Roa-Ureta
AZTI - Tecnalia / Marine Research Unit
Txatxarramendi Ugartea z/g
48395 Sukarrieta (Bizkaia)
SPAIN
Daniel, si nos cuentas brevemente la finalidad de tu análisis exploratorio, podríamos orientarte mejor. ¿Quieres estimar un cuantil? ¿Detectar sub-poblaciones? etc, .... Un saludo. Olivier -- ____________________________________ Olivier G. Nuñez Email: onunez en iberstat.es Tel : +34 663 03 69 09 Web: http://www.iberstat.es ____________________________________ El 09/02/2010, a las 12:01, daniel pacheco gomez escribió:> > 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?. > > He probado hacer varias cosas, aunque no me termina de convencer la > manera de hacerla, ni los datos obtenidos. > > 1 - Una es con la libreria "nortest", que la encontre en un > documento que se llama "FITTING DISTRIBUTIONS WITH R" Ricci, en > ella viene una función que se llama "ad.test()", donde al pasarle > la muestra te devuelve un AD(anderson darling) y p-value de la > distribución, entonces comparo el p-value con la normal, y si > cuadra, pues teoricamente deberia de probar con las otras > distribuciones. Es decir, deberia de calcular los parametros para > la "rweibull" por ejemplo, y aplicarla por ejemplo a 100 > iteraciones, y volver a pasarle el ad.test(), y ver si el p-value > obtenido coincide. Aunque esta forma de hacerlo no me convence, > pues tendria que ir probando una a una cada distribucion, ir > calculandole los parametros y tal, seria una forma muy manual. > Tampoco se si para calcular la distribucion que es, deberia > aplibarle por ejemplo en la weibul, la rweibul, pweibul.. > > 2 - Otra forma de hacerlo, o mas o menos un camino seguir, fue con > la libreria "qAnalyst", la cual tiene un par de funciones para > esto, por ejemplo, la funcion "funInfoFun", donde se le va pasando > la muestra de datos, y la distribucion, y creo, que directamente te > da los valores de AD y p-value, lo malo es que esta funcion solo es > para distribuciones normal, lognormal, gamma, > weibull, logistic, lo cual si no son ninguna de estas no sabria que > hacer, y que al comparar el resultado con el primer metodo no me > coinciden los datos. > > Un ejemplo seria: > > data(warpTiles) > x=warpTiles$warping > infoX=funInfoFun(x,"weibull") > andersonDarlingFun(x=x, infoX$densfun, infoX$theta) > > Lo bueno de esta libreria es para calcular el process capability > segun la distribucion, es decir, las funciones "capabilityNormal" y > "capabilityNotNormal", la cual esta ultima habria que pasarle el > tipo de distribución. > > > > Bueno, como veis, me lio bastante con esto, perdonad por el tocho > escrito, pero era para ser lo mas detallado y descriptivo posible. > > Si alguien ha hecho esto alguna vez, sabe alguna otra manera de > hacerlo, sabe como hacerlo con alguno de estos dos metodos pero > bien, o cualquier cosa, pues le estaria muy agradecio. > > > Asi que gracias de antemano a tod en s. > > > > > > > > > > Navega con el navegador más seguro de todos. ¡Descárgatelo ya! > > _________________________________________________________________ > [[elided Hotmail spam]] > inero! > > [[alternative HTML version deleted]] > > _______________________________________________ > R-help-es mailing list > R-help-es en r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es
> No es Rubén quien pidió los scripts, sino Daniel. > > Ruben cuestionó hacerlo automaticamente como pide Daniel.Si! mis disculpas por la confusion. Pablo> > Un saludo. > > Pablo Emilio Verde escribió: >> 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 >> >