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

Question

我一直在努力将分布与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]"optim中的错误(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函数，但我似乎找不到可以用来实现代码目标的简单示例。有人能帮忙吗？

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

你不会拒绝Normal的零假设。

参考：https://web.archive.org/web/20150407031710/http://worldofpiggy.com:80/2014/02/25/automatic-分布式设置r/

Answer 2

我认为错误主要是因为您的数据。如错误消息中所示，创建NaN使得函数似乎无法获得分数(通过对密度函数进行微分)。[密度函数的范围是非负的，不是吗？]

使用了更简单的矩量法，而不是最大似然估计，它在发出警告的情况下产生参数估计。

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