¿Cómo se hace el bootstrap de la distribución que mejor se ajusta a una muestra?

Question

Si tengo una muestra:

set.seed(0)
x <- rlnorm(500)

Entonces puedo utilizar la función fit.distr para encontrar el mejor ajuste entre dos distribuciones candidatas, por ejemplo

library(MASS)
find.bestfit <- function(x){
   logN <- fitdistr(x, "lognormal")
   gam  <- fitdistr(x, "gamma")
   ans <- ifelse(AIC(logN) < AIC(gam), "logN", "gam")
   return(ans)
}

find.bestfit(x)
[1] "logN"

Sin embargo, existe cierta probabilidad de que no recupere la "verdadera" distribución que se muestreó (en este caso se utilizó "lognormal" para simular x ). ¿Cómo puedo calcular esta probabilidad?

Sólo he llegado a considerar el uso de un enfoque bootstrap, pero no estoy familiarizado con esta técnica y no estoy seguro de por dónde empezar:

## create an empty vector
fit.samps <- rep(NA, 100)
## determine fit to subsamples from original distribution
for(i in 1:100){
  fit.samps[i] <- find.bestfit(sample(x, 10))
}

Sospecho que mi planteamiento es erróneo, porque el tamaño de la muestra es arbitrario, y en última instancia, la distribución mejor ajustada basada en la fitdistr debería seleccionarse la mayoría de las veces.

Agradecería algunas indicaciones sobre cómo podría aplicar el enfoque bootstrap para responder a esta pregunta.

Answer 1

Dado que usted sabe que $X$ es lognormal o gamma, puede utilizar un bootstrap paramétrico en lugar de la versión no paramétrica que usted propuso. Entonces, usted remuestrear a partir de la distribución ajustada y calcular la probabilidad de que find.bestfit da la respuesta correcta.

Esta probabilidad dependerá de si $X$ es lognormal o gamma por lo que hay que hacer dos cálculos distintos.

Aquí hay una manera de hacer esto en R:

library(MASS)

n<-500 # Sample size
B<-100 # Number of bootstrap samples

set.seed(0)
x <- rlnorm(500)

## Create an empty vector
fit.samps <- rep(NA, B)

####

# LOGNORMAL DISTRIBUTION

# Lognormal parameters:
lnpar<-fitdistr(x, "lognormal")$estimate

# Determine fit to parametric bootstrap samples from original distribution
for(i in 1:B){
  fit.samps[i] <- find.bestfit(rlnorm(n,as.vector(lnpar)))
}

# Probability of correct classification if lognormal:
sum(fit.samps=="logN")/B

####

# GAMMA DISTRIBUTION

# Gamma parameters:
gammapar<-fitdistr(x, "gamma")$estimate

##  Determine fit to parametric bootstrap samples from original distribution
for(i in 1:B){
  fit.samps[i] <- find.bestfit(rgamma(n,as.vector(gammapar)))
}

# Probability of correct classification if gamma:
sum(fit.samps=="gam")/B

Para $n=500$ estas probabilidades son ambas prácticamente 1. Para $n\approx 50$ (o menos), se obtienen diferentes probabilidades sin embargo.

Answer 2

El bootstrap se puede utilizar para esto, aunque no se suele hacer. El enfoque sería muestrear con reemplazo n veces de su muestra de tamaño n. Cada vez que se muestre con reemplazo se calculan las estadísticas de bondad de ajuste para las distribuciones que compiten y se elige la distribución que mejor se ajusta. Tome el número de veces que se selecciona la distribución A dividido por el número total de muestras bootstrap para obtener una estimación de la probabilidad de que se seleccione la distribución A.