在R中将伽玛混合物拟合为生育计划 [英] Fit gamma mixture to fertility schedule in R

问题描述

我正在尝试将伽玛混合模型(两个伽玛分布)拟合为年龄-生育能力曲线.我有一个包含特定年龄的生育率和年龄的数据集，我想拟合两个伽玛值以找到相应的参数(最后，我将使用不同年份的生育率概况，并尝试查看这些参数随时间的变化).到目前为止，我尝试使用mixtools库(gammamixEM)，但没有成功.我将非常感谢您的帮助.

I am trying to fit a gamma mixture model (two gamma distributions) to an age-fertility profile. I have a dataset containing age specific fertility rates and age, and I want to fit two gammas in order to find the corresponding parameters (in the end I will use fertility profiles from different years and try to see how the parameters evolve over time). I have so far tried to use mixtools library (gammamixEM) but without success. I would be very grateful for some help.

强麦

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)

推荐答案

之所以无法运行，是因为您的数据集中有零.这是您可以做的:

The reason why it won't run is because you have zeroes in your data set. Here is what you could do:

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

现在您可以适合伽玛混合物

Now you can fit the gamma mixture

require(mixtools)
g3 <- gammamixEM(aa)

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

Now check that it looks OK by plotting the fitted mixture density.

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的伽玛分布参数化为R.为什么呢?谁知道?

Here is another pitfall: gammamixEM apparently parametrises the gamma distribution differently to R. Why? Who knows?

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

看起来很合理，甚至远非完美.

Looks reasonable, if far from perfect.

这篇关于在R中将伽玛混合物拟合为生育计划的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！