[R-es] Una pregunta de estadística (marginalmente relacionada con R)

[Carlos J. Gil Bellosta]

Muchísimas gracias... voy a echarle un vistazo a ver qué tal se comporta...

El día 30 de abril de 2009 18:20, Jorge Ivan Velez
<jorgeivanvelez en gmail.com> escribió:
>
> Carlos,
> Esta es mi propuesta:
> # Librerias
> require(MASS)
> require(survival)
> require(sm)
> # Distribuciones comunes en MASS y survival
> dis <- c("exponential", "lognormal", "logistic", "weibull")
> # Vector de datos
> set.seed(123)
> x <- rweibull(100, shape = 2, scale = 10)
> n <- length(x)
> # Ejemplo para la weibull
> param <- fitdistr(x,dis[4])[[1]]
> simu <- rweibull(n, param[1], param[2])
> # AIC
> tmpf <- extractAIC(survreg(Surv(x)~1, dist=dis[4]))[2]
> tmpf #[1] 579.6235
> # Densidades
> ddx <- density(x)
> ddsimu <- density(simu)
> plot(ddx, xlim=range(ddx$x, ddsimu$x), ylim=range(ddx$y, ddsimu$y))
> points(ddsimu, type='l', col=2)
> legend('topleft',c('Real','Simulados'), col=1:2, lty=1)
>
> # Ahora comparemos las densidades via sm.density.compare en sm
> # Preparando los datos
> y <- c(x, simu)
> gr <- rep(c(1,2), each = n)
> sm.density.compare(y, gr, model="equal")  # Test of equal densities:
>  p-value =  0.56
> Lo que quedaria faltando es la implementación del bootstrap para el AIC,
> pero creo que ya es un poco más fácil. Podrías contruir todo lo anterior en
> una función y luego usar replicate().
> Otra opción es usar ks.test() para comparar la distribución de x con la
> distribución de los datos generados via simulación:
>
> # Esto es SOLO para la weibull
> # -- habria que modificarlo para las demas
> estimeKS <- function(x, alpha=0.05){
>
>       # Estimación de la distribución (dist) y cálculo del AIC
>       n <- length(x)
>       param <- fitdistr(x,dis[4])[[1]]
>       simu <- rweibull(n, param[1], param[2])
>       # KS
>       tmpf <- ks.test(x, simu)$p.value
>       tmpf > alpha  # Acepto Ho: F(x) = G(simu)?
> }
> # --------
> # Ejemplo
> # --------
> # Datos
> set.seed(123)
> x <- rweibull(100, shape = 2, scale = 10)
> KSs <- replicate(1000, estimeKS(x))
> sum(KSs)/1000   # [1] 0.998
>
> Un saludo,
> Jorge Ivan Velez
>
>
> 2009/4/30 Jorge Ivan Velez <jorgeivanvelez en gmail.com>
>>
>> Hola Carlos,
>> Podrías hacer las simulaciones con diferentes tamaños de muestra (n)
>> cuando generas las observaciones de la distribución específica. Lo común es
>> usar n como la longitud del vector que quieres ajustar a determinada
>> distribución; sin embargo, esto no quiere decir que no puedes utilizar algún
>> otro valor.
>> Ten en cuenta que bootstrap es importante para validar que efectivamente
>> la distribución que ajustas es "la que potencialmente" generó los datos como
>> dice Pablo.
>>
>> En cuanto al programa, sería algo como lo siguiente ( en palabras :-(  ):
>> 1. Tomar el vector de observaciones (digamos x) y ajustar una distribución
>> (digamos F*) dentro de la gama de posibilidades usando, por ejemplo,
>> fitdistr en MASS.
>> 2. Generar n (o más) observaciones de F* y calcular el AIC. (Deberías usar
>> el código de Pablo :-) ).
>> 3. Repetir el paso anterior N veces y construir, por ejemplo, un intervalo
>> de confianza para el AIC estimado. Para ello puedes usar la función boot en
>> la libreria boot.
>> 4. Repetir 1-3 para las demás distribuciones en tu gamma.
>> 5. Reportar los resultados y comparar.
>> Pablo: Muchas gracias por el código que envías; es fácil de implementar y
>> entender. Alguna idea de cómo extenderlo, por ejemplo a la distribución
>> gamma?
>> Un saludo a ambos,
>> Jorge Ivan Velez
>>
>> 2009/4/30 Carlos J. Gil Bellosta <cgb en datanalytics.com>
>>>
>>> Muchas gracias por la contestación (y bienvenido a la lista).
>>>
>>> Pensé en utilizar AIC pero me da algo de miedo cuando los modelos no
>>> están anidados (la exponencial lo está, en cierto modo, en la
>>> weibull), pero, en general...
>>>
>>> ¿Has probado a experimentar con tu código y a probar distintos valores
>>> del número de datos y de distribuciones con los que los generas?
>>>
>>> Un saludo,
>>>
>>> Carlos J. Gil Bellosta
>>> http://www.datanalytics.com
>>>
>>> El día 30 de abril de 2009 16:23, Pablo Emilio Verde
>>> <PabloEmilio.Verde en uni-duesseldorf.de> escribió:
>>> > Este ejemplo te puede servir. Utilizo las distribuciones que estan la
>>> > funcion survreg() del paquete survival y extraigo en AIC con la
>>> > funcion extractAIC() del paquete MASS.
>>> >
>>> > #################################################
>>> > set.seed(123)
>>> > x <- rweibull(100, shape = 2, scale = 10)
>>> >
>>> > library(MASS) # para aplicar extractAIC
>>> > library(survival) # para survreg
>>> >
>>> > distrib <-c("weibull", "exponential", "gaussian", "logistic",
>>> >  "lognormal", "loglogistic")
>>> >
>>> > for( Dis in distrib){ tmpf <-extractAIC(survreg(Surv(x)~1, dist=Dis))
>>> >                             cat(Dis, ", AIC = ",tmpf[2], "\n")
>>> >               }
>>> > #################################################
>>> >
>>> > El menor valor de AIC indica la distribucion de probabiliad que
>>> > potencialmente genero los datos.
>>> >
>>> > Pablo
>>>
>>> _______________________________________________
>>> R-help-es mailing list
>>> R-help-es en r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>
>
>