r - 如何在 R 中将分布拟合到样本数据？

Question

我一直在努力将分布拟合到 R 中的样本数据。我已经考虑过使用 fitdist 和 fitdistr 函数，但我似乎在使用这两个函数时都遇到了问题。

快速背景；我的代码的输出应该是带有参数的最适合所提供数据的分布(来自分布列表)。这需要在没有人类交互的情况下发生，因此比较图表不是一种选择。我想我可以将每个分布拟合到数据，从卡方检验中得出 p 值，并找到具有最高 p 值的分布。我在样本数据的正态分布方面取得了一些成功，但是一旦我尝试拟合更复杂的东西( Gamma 分布，如代码中所示)，我就会遇到各种错误。我究竟做错了什么？

library(fitdistrplus) 
require(MASS) 
set.seed(1) 
testData <- rnorm(1000) 
distlist <- c("norm","unif","exp")

(z <- fitdist(testData,"gamma",start=list(rate=0.1),fix.arg=list(shape=4)))

我遇到的错误示例是:

[1] "Error in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data, : \n initial value in 'vmmin' is not finite\n" attr(,"class")

和

Error in fitdist(testData, "gamma", start = list(rate = 0.1), fix.arg = list(shape = 4)) : the function mle failed to estimate the parameters, with the error code 100

我知道我可能错误地实现了 fitdist 函数，但我似乎找不到可以适应以实现我的代码目标的简单示例。有人可以帮忙吗？

Answer 1

您正在寻找 Kolmogorov-Smirnov 检验。零假设是指数据样本来自假设的分布。

fitData <- function(data, fit="gamma", sample=0.5){
 distrib = list()
 numfit <- length(fit)
 results = matrix(0, ncol=5, nrow=numfit)

 for(i in 1:numfit){
if((fit[i] == "gamma") | 
     (fit[i] == "poisson") | 
     (fit[i] == "weibull") | 
     (fit[i] == "exponential") |
     (fit[i] == "logistic") |
     (fit[i] == "normal") | 
     (fit[i] == "geometric")
) 
  distrib[[i]] = fit[i]
else stop("Provide a valid distribution to fit data" )
 }

 # take a sample of dataset
 n = round(length(data)*sample)
 data = sample(data, size=n, replace=F)

 for(i in 1:numfit) {
  if(distrib[[i]] == "gamma") {
  gf_shape = "gamma"
  fd_g <- fitdistr(data, "gamma")
  est_shape = fd_g$estimate[[1]]
  est_rate = fd_g$estimate[[2]]

  ks = ks.test(data, "pgamma", shape=est_shape, rate=est_rate)

  # add to results
  results[i,] = c(gf_shape, est_shape, est_rate, ks$statistic, ks$p.value)
}

else if(distrib[[i]] == "poisson"){
  gf_shape = "poisson"
  fd_p <- fitdistr(data, "poisson")
  est_lambda = fd_p$estimate[[1]]

  ks = ks.test(data, "ppois", lambda=est_lambda)
  # add to results
  results[i,] = c(gf_shape, est_lambda, "NA", ks$statistic, ks$p.value)
}

else if(distrib[[i]] == "weibull"){
  gf_shape = "weibull"
  fd_w <- fitdistr(data,densfun=dweibull,start=list(scale=1,shape=2))
  est_shape = fd_w$estimate[[1]]
  est_scale = fd_w$estimate[[2]]

  ks = ks.test(data, "pweibull", shape=est_shape, scale=est_scale)
  # add to results
  results[i,] = c(gf_shape, est_shape, est_scale, ks$statistic, ks$p.value) 
}

else if(distrib[[i]] == "normal"){
  gf_shape = "normal"
  fd_n <- fitdistr(data, "normal")
  est_mean = fd_n$estimate[[1]]
  est_sd = fd_n$estimate[[2]]

  ks = ks.test(data, "pnorm", mean=est_mean, sd=est_sd)
  # add to results
  results[i,] = c(gf_shape, est_mean, est_sd, ks$statistic, ks$p.value)
}

else if(distrib[[i]] == "exponential"){
  gf_shape = "exponential"
  fd_e <- fitdistr(data, "exponential")
  est_rate = fd_e$estimate[[1]]
  ks = ks.test(data, "pexp", rate=est_rate)
  # add to results
  results[i,] = c(gf_shape, est_rate, "NA", ks$statistic, ks$p.value)
}

else if(distrib[[i]] == "logistic"){
  gf_shape = "logistic"
  fd_l <- fitdistr(data, "logistic")
  est_location = fd_l$estimate[[1]]
  est_scale = fd_l$estimate[[2]]
  ks = ks.test(data, "plogis", location=est_location, scale=est_scale)
  # add to results
  results[i,] = c(gf_shape, est_location, est_scale, ks$statistic,    ks$p.value) 
    }
  }
  results = rbind(c("distribution", "param1", "param2", "ks stat", "ks    pvalue"),   results)
  #print(results)
  return(results)
  }

应用于您的示例:

library(MASS)
set.seed(1) 
testData <- rnorm(1000) 
res = fitData(testData, fit=c("logistic","normal","exponential","poisson"),
    sample=1)
res

您不拒绝正态的原假设。

引用:https://web.archive.org/web/20150407031710/http://worldofpiggy.com:80/2014/02/25/automatic-distribution-fitting-r/