使用 mle2 函数

Question

我想在这样的模型中找到参数 epsilon 和 mu 的 MLE：

$$X \sim \frac{1}{mu1}e^{-x/mu1}+\frac{1}{mu2}e^[-x/mu2}$$

library(Renext)
library(bbmle)
epsilon = 0.01

#the real model
X <- rmixexp2(n = 20, prob1 = epsilon, rate1 = 1/mu1, rate2 = 1/mu2)

LL <- function(mu1,mu2, eps){
  R = (1-eps)*dexp(X,rate=1/mu1,log=TRUE)+eps*dexp(X,rate=1/mu2,log=TRUE)
  -sum(R)
}
fit_norm <- mle2(LL, start = list(eps = 0,mu1=1, mu2 = 1), lower = c(-Inf, 0),
                 upper = c(Inf, Inf), method = 'L-BFGS-B')

summary(fit_norm)


But I get the error 

> fn = function (p) ':method 'L-BFGS-B' requires finite values of fn"

Answer 1

这里有很多问题。第一个是您的可能性表达式是错误的（您不能单独记录组件然后添加它们，您必须添加组件并然后记录日志）。你的界限也很有趣：混合概率应该是 [0,1]，平均值应该是 [0, Inf]。

您遇到的另一个问题是，对于当前的模拟设计 (n=20, prob=0.01)，您很有可能在第一个混合组件中获得 no 个点（一个点在第二个分量中是 1-0.01=0.99，所以 所有 个点在第二个分量中的概率是 0.99^20 = 82%）。在这种情况下，MLE 将是退化的（即，您试图将双组分混合物拟合到基本上只有一个组分的数据集）；在这种情况下，这些解决方案中的任何一个都会给出等效的可能性：

• prob=0，mu2=数据的平均值，mu1=anything
• prob=1，mu1=数据的平均值，mu2=anything
• mu1=mu2=数据的平均值，prob=anything

使用所有这些解决方案，最终的结果将非常敏感地取决于起始条件和优化算法。

对于这个问题，我鼓励您使用Renext 包中的内置 dmixexp2 函数（它正确地将对数似然实现为log(p*Prob(X|exp1) + (1-p)*Prob(X|exp2))）和公式接口给mle2：

fit_norm <- mle2(X ~ dmixexp2(rate1=1/mu1,rate2=1/mu2,prob1=eps),
                 data=list(X=X),
                 start = list(mu1=1, mu2 = 2, eps=0.4),
                 lower = c(mu1=0, mu2=0, eps=0),
                 upper = c(mu1=Inf, mu2=Inf, eps=1),
                 method = 'L-BFGS-B')

这给了我mu1=1.58、mu2=2.702、eps=0 的估计值。 mean(X) 在我的例子中等于 mu2 的值，所以这是上面项目符号列表中的第一个例子。您还会收到警告：

某些参数在边界上：基于 Hessian 的方差-协方差计算可能不可靠

还有多种更专业的算法用于拟合混合模型（尤其是那些基于期望最大化算法的算法）；您可以在 CRAN 上查找包（flexmix 就是其中之一）。

这个问题足够小，您可以通过蛮力（下面的代码）可视化整个对数似然面：颜色表示与最小负对数似然的偏差（颜色渐变是对数缩放的，所以有一个小的偏移以避免 log(0))。深蓝色代表最适合数据的参数，黄色是最差的。

---

dd <- expand.grid(mu1=seq(0.1,4,length=51),
                  mu2=seq(0.1,4,length=51),
                  eps=seq(0,1,length=9),
                  nll=NA)
for (i in 1:nrow(dd)) {
    dd$nll[i] <- with(dd[i,],
                      -sum(dmixexp2(X,rate1=1/mu1,
                                    rate2=1/mu2,
                                    prob1=eps,
                                    log=TRUE)))
}

library(ggplot2)
ggplot(dd,aes(mu1,mu2,fill=nll-min(nll)+1e-4)) +
    facet_wrap(~eps, labeller=label_both) +
    geom_raster() +
    scale_fill_viridis_c(trans="log10") +
    scale_x_continuous(expand=c(0,0)) +
    scale_y_continuous(expand=c(0,0)) +
    theme(panel.spacing=grid::unit(0.1,"lines"))
ggsave("fit_norm.png", type="cairo-png")