Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

If I have a sample:

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

Then I can use the fit.distr function to find the best fit among two candidate distributions, e.g.

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"

However, there is some probability that I will not recover the "true" distribution that was sampled (in this case "lognormal" was used to simulate x). How can I calculate this probability?

I have only gotten so far as to consider using a bootstrap approach, but I am not familiar with this technique and am not sure exactly where to start:

## 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))
}

I suspect that my approach is wrong, because the sample size is arbitrary, and ultimately, the best fit distribution based on the fitdistr function should be selected most of the time.

I would appreciate some pointers on how I might apply the bootstrap approach to answer this question.

share|improve this question

2 Answers 2

up vote 5 down vote accepted

Since you know that $X$ is either lognormal or gamma, you can use a parametric bootstrap instead of the nonparametric version that you proposed. You would then resample from the fitted distribution instead, and compute the probability that find.bestfit gives the right answer.

This probability will depend on whether $X$ is lognormal or gamma, so you have to make two separate computations.

Here is a way to do this in 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

For $n=500$ these probabilities are both virtually 1. For $n\approx 50$ (or less), you get different probabilities though.

share|improve this answer

The bootstrap can be used for this although it is not commonly done. Th approach would be to sample with replacement n times from your sample of size n. Each time you sample with replacement you compute the goodness of fit statistics for the competing distributions and pick the distribution that fits best. Take the number of times distribution A is selected divided by the total number of bootstrap samples to get an estimate for the probability that distribution A will be selected.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.