如何将分布拟合到 R 中的样本数据？

Posted 2023-02-19

【中文标题】如何将分布拟合到 R 中的样本数据？

【英文标题】：How do I fit distributions to sample data in R?

【发布时间】：2015-07-03 13:23:39

【问题描述】：

我一直在努力将分布拟合到我在 R 中拥有的样本数据。我已经研究过使用 fitdist 和 fitdistr 函数，但我似乎都遇到了问题。

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

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] "优化错误(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data, : \n 'vmmin' 中的初始值不是有限的\n" attr(,"class")

和

fitdist 中的错误(testData, "gamma", start = list(rate = 0.1), fix.arg = list(shape = 4)) ：函数 mle 无法估计参数， 错误代码 100

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

【问题讨论】：

错误消息说明了一切：对数似然在初始值处不是有限的。伽马分布有正支持，而样本肯定有负值，因此对数似然是无限的。

嗯。甚至从未考虑过这一点；你说的对。我将尝试对样本数据进行一些控制，以仅包含正数据。感谢您的反馈，伙计。

密切相关：stats.stackexchange.com/questions/30491/…，***.com/questions/2661402/…

另外，我不建议使用 p 值进行模型选择，它们不表示观察结果由特定模型生成的概率。 Akaike information criterion 将是一个简单、易于计算的替代方案。

@Arpi，非常感谢您的建议。我将阅读该技术，看看它是否效果更好。任何帮助或建议都受到高度重视，因此我非常感谢。

【参考方案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/

【讨论】：

哇哦。非常感谢你，它准确地解释了我需要什么。我假设我可以扩展此功能以包括其他发行版？抱歉，如果我的问题听起来很愚蠢，我对 R 的了解仍在增长。

是的，您可以将其扩展到其他发行版，但您对数据性质的了解应该会将您的选择限制在几个发行版中。

KS 检验不适用于拟合分布，即从数据估计分布的参数。在这种情况下，KS 统计数据不是无房客分布的，并且临界值不准确。有一些技术可以修改统计数据，即使在估计参数的情况下，它也确实变得无分布，但它不是微不足道的（例如，参见Khmaladze Transform）。

带有估计参数的 Kolmovorov-Smirnov-Test (KS-Test) 的 p 值将是非常错误的。所以不幸的是，您不能只拟合一个分布，然后使用 Kolmogorov-Smirnov-Test 中的估计参数来测试您的样本。看这个帖子stats.stackexchange.com/questions/132652/…

【参考方案2】：

我认为错误主要是因为您的数据。如错误消息中所示，创建了NaN，因此该函数似乎无法获得分数（通过区分密度函数）。 [密度函数的范围是非负的，不是吗？]

Method of moments 更简单，用于代替最大似然估计，尽管有警告，它仍会生成参数估计。

library(fitdistrplus) 
require(MASS) 
set.seed(1) 
testData <- rnorm(1000) 
fitdist(testData, "gamma", method = "mme", start = list(shape = 0.1, rate = 0.1))

Fitting of the distribution ' gamma ' by matching moments 
Parameters:
           estimate
shape  0.0001268054
rate  -0.0108863200
Warning message:
In dgamma(c(-0.626453810742332, 0.183643324222082, -0.835628612410047,  :
  NaNs produced

【讨论】：

非常感谢您的详细回答。实际上，我从未考虑过使用 MME，正如您所指出的，它非常适合这个问题。我假设我可以使用 gofstat 之类的方法简单地从中得出卡方 p 值？

MME 仅使用矩来拟合分布，而 MLE 通过拟合似然函数使用更多信息，我想这就是为什么前者至少返回结果的原因。 |伽马分布的域是 [0, infinity) 而对于正态分布它是 (-infinity, infinity)，因此随机样本的负实现会导致问题。因此，如果严格控制值，我猜fitdist() 应该返回一个结果。 |第二种方式也可以作为非参数估计的一种方式。