classFitDist <- c(class(fe), class(flogis), class(fn), class(fp),class(fg))

distributions <- classFitDist == "fitdist"


  AIC <- data.frame()
  if(class(fe)=="fitdist") {AIC[1,ncol(AIC)+1] <- fe$aic}
  if(class(flogis)=="fitdist") {AIC[1,ncol(AIC)+1] <- flogis$aic}
  if(class(fn)=="fitdist") {AIC[1,ncol(AIC)+1] <- fn$aic}
  if(class(fp)=="fitdist") {AIC[1,ncol(AIC)+1] <- fp$aic}
  if(class(fg)=="fitdist") {AIC[1,ncol(AIC)+1] <- fg$aic}

  names(AIC) <- c("exp", "logis", "norm", "pois", "gamma")[distributions]

  fit <- names(AIC[which.min(AIC)])


  mean <- switch (fit,
                  exp = 1/fe$estimate[[1]],
                  logis = flogis$estimate[[1]],
                  norm = fn$estimate[[1]],
                  pois = fp$estimate[[1]],
                  gamma = fg$estimate[[1]]/fg$estimate[[2]]
  )

  sd <- switch (fit,
                exp = mean,
                logis = (flogis$estimate[[2]]*pi)/sqrt(3),
                norm = fn$estimate[[2]],
                pois = sqrt(mean),
                gamma = sqrt(fg$estimate[[1]]/(fg$estimate[[2]]^2))
  )

  # etc, etc.

goodnessoffit <- gofstat(list(fe, flogis, fn, fp, fg)[distributions],  fitnames = c("exp", "logis", "norm", "pois","gamma")[distributions])
fit <- names(which(goodnessoffit$aic == min(goodnessoffit$aic)))