将伽马混合物拟合到 R 中的生育计划答案

【问题标题】：Fit gamma mixture to fertility schedule in R将伽马混合物拟合到 R 中的生育计划
【发布时间】：2013-07-01 07:07:02
【问题描述】：

我正在尝试将伽马混合模型（两个伽马分布）拟合到年龄-生育力曲线。我有一个包含特定年龄生育率和年龄的数据集，我想拟合两个伽马以找到相应的参数（最后我将使用不同年份的生育率资料并尝试查看参数如何随时间演变）。到目前为止，我尝试使用 mixtools 库（gammamixEM）但没有成功。我将非常感谢您的帮助。

麦芽酒

a<- structure(list(EDAD = structure(1:45, .Label = c("11", "12", "13", "14", "15",   
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", 
"30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", 
"44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "Total" ),
class = "factor"), value = c(0, 0, 0, 0, 0.002761668, 0.006712018, 0.010820244, 
0.017867778, 0.029533765, 0.034055242, 0.036665669, 0.043840421, 0.042949584, 
0.042344612, 0.050588917, 0.050187588, 0.054114728, 0.057258792, 0.059280324, 
0.062566731, 0.062369629, 0.062154767, 0.063734337, 0.058236776, 0.052623842, 
0.046330921, 0.040639027, 0.033707865, 0.02531141, 0.017651534, 0.010953808, 
0.007463863, 0.003224766, 0.002190101, 0.001117443, 0.000465116, 0.000363901, 
0.00012647, 0.000267326, 0.000280308, 0, 0, 0, 0, 0)), .Names = c("EDAD", "value"), 
class = "data.frame", row.names = 79596:79640)

【问题讨论】：

你有什么问题？
如何使伽玛混合物适应我的生育计划

标签： r gamma-distribution

【解决方案1】：

它无法运行的原因是因为您的数据集中有零。您可以这样做：

aa <- a$value[a$value > 0]

现在您可以拟合伽马混合物

require(mixtools)
g3 <- gammamixEM(aa)

现在通过绘制拟合的混合密度来检查它看起来是否正常。

d3 <- function(x) g3$lambda[1]*dgamma(x, g3$gamma.pars[1], 1/g3$gamma.pars[2]) + g3$lambda[2]*dgamma(x, g3$gamma.pars[3], 1/g3$gamma.pars[4])

这是另一个陷阱：gammamixEM 显然对伽马分布的参数化与 R 不同。为什么？谁知道？

x <- seq(min(aa), max(aa), 0.001)
plot(x, d3(x), "l")
hist(aa, col="pink", add=T, freq=F, breaks=10)

看起来很合理，但远非完美。

【讨论】：