如何在R中将伽玛分布拟合到数据中？

Question

假设我有一个变量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

我该如何将伽玛分布适配到它？

- user1172468

---

请参阅函数MASS :: fitdistr。 - Rui Barradas

抱歉，我刚才更仔细地看了你的数据和它的图形。你为什么想要对其进行伽马分布拟合？all.equal(seq(0, 1, by = 0.01), x)的返回值是TRUE。 - Rui Barradas

@RuiBarradas -- 谢谢 -- 让我更新问题 -- 我通过向伽马偏差添加少量高斯噪声来生成数据 -- 让我在注释中发布这个。 - user1172468

@RuiBarradas -- 感谢您的更新问题。 - user1172468

在我的书（第433页）中有一个对于高斯污染伽马样本进行MLE的例子，尽管我怀疑它不适用于这个数据集——高斯分量的标准差可能太小而无法区分它。 - Ben Bolker

Answer 1

一个很好的选择是由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)

- Edgar Santos

---

这个解决方案可能存在问题。x有负值。sum(x < 0)返回5。将MASS::fitdistr应用于此向量时会出现错误：Error in fitdistr(x, "gamma") : gamma values must be >= 0。因此，也许fitdistrplus包正在执行它应该执行的所有检查。 - Rui Barradas

@RuiBarradas，谢谢。很好的建议。已经进行了编辑。该函数不接受负值。 - Edgar Santos

Answer 2

您可以尝试快速适配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
> 

- Severin Pappadeux

---

这种方法被称为“矩法”，但我希望在修补裂缝时非常小心（因为所谓的伽马分布样本中存在负值）。 - Ben Bolker

Answer 3

您还可以尝试使用OneStep包中的onestep命令，使用快速高效的Le Cam单步估计过程来拟合Gamma分布。

library(OneStep)

set.seed(2017)
x <- rgamma(100,2,11) + rnorm(100,0,.01)
onestep(x,"gamma")

Parameters:
       estimate
shape  2.183281
rate  12.837968

- Alexandre Brouste