1
$\begingroup$

I'm trying to write a piece of code in R that

  • finds the most-fitting distribution to a set of data, by
  • performing goodness-of-fit tests to a list of distributions, and then
  • finding the most fitting one.
  • This program should be able to run in real-time, so analysis should be very light on computational load. What I mean by this is that it should be able to process, say, a fit every second or few seconds at the most, so the simpler the program, the better.

For instance, I've already written the following code:

for(i in 1:numfit) {
if(distrib[[i]] == "negative binomial"){
  gf_shape = "negative binomial"
  fd_nb <- tryCatch((fitdistr(data, "negative binomial", start=list(size=1,prob=0.5))),
    error = function(fd_nb) fd_nb <- fitdistr(data, "negative binomial"))
  est_size = fd_nb$estimate[[1]]
      est_prob = fd_nb$estimate[[2]]
  gfn = goodfit(data,type="nbinomial",method="MinChisq",par = list(size = est_size))
  tidied = tidy(summary(gfn))
  results[i,] = c(gf_shape, est_lambda, "NA", tidied$X.2, tidied$P...X.2.)
}

else if(distrib[[i]] == "poisson"){
  gf_shape = "poisson"
  fd_p <- fitdistr(data, "poisson")
  est_lambda = fd_p$estimate[[1]]
  gf = goodfit(data,type="poisson",method="MinChisq",par = list(lambda = est_lambda))
  tidied = tidy(summary(gf))
  results[i,] = c(gf_shape, est_lambda, "NA", tidied[1,1], tidied[1,3])
}
results = rbind(c("distribution", "parameter 1", "parameter 2", "chi-squared test statistic", "P > X2"),   results)
return(results)
}

that performs a chi-square goodness-of-fit test of my data to a Poisson and negative binomial distribution, from which I can then find the distribution with the lowest chi-squared test statistic and infer the most suitable distribution to the data from there.

My question is how to do this with continuous data/continuous distributions. Using the goodfit package in R worked really well for me, but it only works with discrete distributions. I do, however, need to use the chi-square goodness-of-fit method (project requirement), so I'm not sure which package or method to turn to and how to implement this.

Can anyone help me with an easy idea of how to implement something similar for, say, a normal distribution? At the moment I think I'm just going to write a piece of code to bin the data into categories, the way you would do the chi-square test by hand, but any way to optimize this process would help. Help would be much appreciated.

$\endgroup$
2
  • $\begingroup$ I should clarify that you'll need to reword your question for the other site, but with some work, it might be a good fit. $\endgroup$ Aug 10, 2015 at 16:26
  • $\begingroup$ I added "continuous" to the title so that this wouldn't sound like a boring introductory statistics question. $\endgroup$
    – EdM
    Aug 10, 2015 at 21:50

2 Answers 2

3
$\begingroup$

It's a pity that a chi-square goodness-of-fit test is a "project requirement," because in general there is much to be lost by binning continuous variables. If possible, try to convince those in charge to allow methods more appropriate for continuous variables in this context of distribution fitting.

If you are stuck with binning you might have a problem. Essentially, if you define the bin boundaries based on parameters of a distribution estimated from the data, then the chi-square test itself may no longer be valid. This is beyond the simple loss of degrees of freedom in chi-square due to estimating parameters; the chi-square statistic as usually calculated may be no longer distributed as chi-square. See the "Level 3" section of the answer by @cardinal on that linked page; that answer includes some recommendations for how to proceed (which are beyond my particular expertise). As you proceed, you do have to consider what you are trying to accomplish by fitting different types of distributions to the same data, and to recognize that results might be hard to generalize even if you find a reliable way to define useful bin boundaries.

$\endgroup$
3
  • $\begingroup$ Awesome, thanks very much for the recommendation. I might be able to convince them otherwise, in which case a KS-test would probably be the easiest way to go, right? $\endgroup$
    – Martin
    Aug 10, 2015 at 22:16
  • $\begingroup$ Easiest perhaps, but you have to take care when using parameter values derived from the data, as explained in the second linked page in my answer. $\endgroup$
    – EdM
    Aug 11, 2015 at 0:12
  • $\begingroup$ For a parallel construction with your discrete-distribution case, you could use the Akaike Information Criterion to compare fits generated by the fitdistr function to continuous distributions. That's analogous to your use of chi-square for discrete distributions, and avoids the difficulties in K-S tests on distributions with parameters estimated from the data.. Again, see my second linked page for examples. I repeat my caution about thinking hard about what you are trying to accomplish in this process. $\endgroup$
    – EdM
    Aug 11, 2015 at 0:24
0
$\begingroup$

One thing is to perform a goodness-of-fit test and another to use it to select an appropriate distribution for the data. For the former, I wrote the following function:

#' Chi-squared Goodness-of-fit Test for Continuous Data
#'
#' Performs the chi-squared goodness-of-fit test for a continuous distribution
#' by grouping data into *bins* (also called *intervals*, *classes* or *cells*)
#' with equal probabilities under the null hypothesis.
#' @param x numeric vector containing the observed values.
#' @param distribution character string naming a continuous distribution,
#' such as "norm" or "unif" (the cumulative distribution `p<distribution>()` and
#' quantile `q<distribution>()` functions must exist).
#' @param nclasses number of bins.
#' @param output logical; if `TRUE` an histogram is plotted and a table with
#' the results for each class is printed.
#' @param nestpar number of estimated parameters (composite null hypothesis).
#' @param ...  parameters of the distribution (specified by `distribution`).
#' @return A list with class `"htest"` containing the following components:
#' \item{statistic}{the value the chi-squared test statistic.}
#' \item{parameter}{the degrees of freedom of the approximate chi-squared
#' distribution of the test statistic.}
#' \item{p.value}{the p-value for the test.}
#' \item{method}{a character string indicating the type of test performed.}
#' \item{data.name}{a character string with the actual `x` argument name.}
#' \item{classes}{a character vector with the class labels.}
#' \item{observed}{the observed counts.}
#' \item{expected}{the expected counts under the null hypothesis.}
#' \item{residuals}{the Pearson residuals, `(observed - expected) / sqrt(expected)`.}
#' @examples
#' nx <- 30
#' x <- rnorm(nx)
#' chisq.cont.test(x, distribution = "norm", nestpar = 2,
#'                 mean = mean(x), sd = sqrt((nx - 1) / nx) * sd(x))
#' @export
chisq.cont.test <- function(x, distribution = "norm",
                            nclasses = floor(length(x)/5), output = TRUE, nestpar = 0, ...) {
  # Funciones distribución
  q.distrib <- eval(parse(text = paste("q", distribution, sep = "")))
  d.distrib <- eval(parse(text = paste("d", distribution, sep = "")))
  # Puntos de corte
  q <- q.distrib((1:(nclasses - 1))/nclasses, ...)
  tol <- sqrt(.Machine$double.eps)
  xbreaks <- c(min(x) - tol, q, max(x) + tol)
  # Gráficos y frecuencias
  if (output) {
xhist <- hist(x, breaks = xbreaks, freq = FALSE,
              lty = 2, border = "grey50")
curve(d.distrib(x, ...), add = TRUE)
  } else {
xhist <- hist(x, breaks = xbreaks, plot = FALSE)
  }
  # Cálculo estadístico y p-valor
  O <- xhist$counts  # Equivalente a table(cut(x, xbreaks)) pero más eficiente
  E <- length(x)/nclasses
  DNAME <- deparse(substitute(x))
  METHOD <- "Pearson's Chi-squared test"
  STATISTIC <- sum((O - E)^2/E)
  names(STATISTIC) <- "X-squared"
  PARAMETER <- nclasses - nestpar - 1
  names(PARAMETER) <- "df"
  PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
  # Preparar resultados
  classes <- format(xbreaks)
  classes <- paste("(", classes[-(nclasses + 1)], ",", classes[-1], "]",
                   sep = "")
  RESULTS <- list(classes = classes, observed = O, expected = E,
                  residuals = (O - E)/sqrt(E))
  if (output) {
    cat("\nPearson's Chi-squared test table\n")
    print(as.data.frame(RESULTS))
  }
  if (any(E < 5))
    warning("Chi-squared approximation may be incorrect")
  structure(c(list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL,
                   method = METHOD, data.name = DNAME), RESULTS), class = "htest")
}

which could be used as a basis for the latter...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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