3

我正在尝试将伽马混合模型（两个伽马分布）拟合到年龄生育率曲线。我有一个包含特定年龄生育率和年龄的数据集，我想拟合两个伽马以找到相应的参数（最后我将使用不同年份的生育率资料并尝试查看参数如何随时间演变）。到目前为止，我尝试使用 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)
4

1 回答 1

3

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

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)

看起来很合理，如果远非完美的话。

在此处输入图像描述

于 2013-07-25T02:52:38.523 回答