您如何将伽马分布拟合到 R 中的数据？答案

【问题标题】：How would you fit a gamma distribution to a data in R?您如何将伽马分布拟合到 R 中的数据？
【发布时间】：2018-01-14 03:06:41
【问题描述】：

假设我有使用以下方法生成的变量 x：

x <- rgamma(100,2,11) + rnorm(100,0,.01) #gamma distr + some gaussian noise

    head(x,20)
 [1] 0.35135058 0.12784251 0.23770365 0.13095612 0.18796901 0.18251968
 [7] 0.20506117 0.25298286 0.11888596 0.07953969 0.09763770 0.28698417
[13] 0.07647302 0.17489578 0.02594517 0.14016041 0.04102864 0.13677059
[19] 0.18963015 0.23626828

我怎样才能将伽马分布拟合到它？

【问题讨论】：

见函数MASS::fitdistr。
抱歉，我刚刚更仔细地查看了您的数据及其图表。你为什么要给它加一个伽马？ all.equal(seq(0, 1, by = 0.01), x) 返回TRUE。
@RuiBarradas -- 谢谢 -- 让我更新一下问题 -- 我通过在伽马偏差中添加少量高斯噪声来生成数据 -- 让我把它贴在笔记中
@RuiBarradas -- 更新了这个问题，谢谢。
在my book（第 433 页）中有一个对这种受高斯污染的 Gamma 样本进行 MLE 的示例，尽管我怀疑它是否适用于该数据集——标准高斯分量的 dev 是概率。太小了，无法区分。

标签： r distribution gamma-distribution

【解决方案1】：

您可以尝试快速拟合 Gamma 分布。作为双参数分布，可以通过找到样本均值和方差来恢复它们。在这里，只要平均值为正，您就可以有一些样本为负。

set.seed(31234)
x <- rgamma(100, 2.0, 11.0) + rnorm(100, 0, .01) #gamma distr + some gaussian noise
#print(x)

m <- mean(x)
v <- var(x)

print(m)
print(v)

scale <- v/m
shape <- m*m/v

print(shape)
print(1.0/scale)

对我来说它打印出来

> print(shape)
[1] 2.066785
> print(1.0/scale)
[1] 11.57765
>

【讨论】：

true（这种方法称为“矩量法”）但我希望非常小心地掩盖裂缝（据称 Gamma 中存在负值-distributed sample) 以这种方式...

【解决方案2】：

ML Delignette-Muller 等人的 fitdistrplus 包是一个不错的选择。例如，使用您的方法生成数据：

set.seed(2017)
x <- rgamma(100,2,11) + rnorm(100,0,.01)
library(fitdistrplus)
fit.gamma <- fitdist(x, distr = "gamma", method = "mle")
summary(fit.gamma)

Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters : 
       estimate Std. Error
shape  2.185415  0.2885935
rate  12.850432  1.9066390
Loglikelihood:  91.41958   AIC:  -178.8392   BIC:  -173.6288 
Correlation matrix:
          shape      rate
shape 1.0000000 0.8900242
rate  0.8900242 1.0000000


plot(fit.gamma)

【讨论】：

此解决方案可能存在问题。 x 有负值。 sum(x < 0)returns 5. MASS::fitdistr 应用于此向量时会出现错误：Error in fitdistr(x, "gamma") : gamma values must be >= 0。所以也许包fitdistrplus 正在做它应该做的所有检查。
@RuiBarradas，谢谢。好点子。上面编辑。该函数不接受负值。