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