R：使用优化的指数混合的最大似然估计答案

【问题标题】：R: Maximum Likelihood Estimation of a exponential mixture using optimR：使用优化的指数混合的最大似然估计
【发布时间】：2020-08-11 06:41:00
【问题描述】：

我正在尝试使用 R 中的对数似然函数和 optim 函数从混合双指数模型中获取参数 w, lambda_1, lambda_2 和 p。模型如下

这里是代码

biexpLL <- function(theta, y) {
  # define parameters
  w <- theta[1]
  lambda_1 <- theta[2]
  a <- theta[3]
  lambda_2 <- theta[4]
  # likelihood function with dexp
  l <- w * dexp((y - a), rate = 1/lambda_1) + (1 - w) * dexp((y - a), rate = 1/lambda_2)
  
  - sum(log(l))
}
# Generate some fake data
w <- 0.7
n <- 500
lambda_1 <- 2
lambda_2 <- 0.2
set.seed(45)
biexp_data <- (w * rexp(n, 1/lambda_1) + (1 - w) * rexp(n, 1/lambda_2)) 
# Optimization
optim(par = c(0.5,0.1,0.001,0.2),
      fn=biexpLL,
      y=biexp_data)
#$par
#[1] -94789220.4     16582.9   -333331.7 134744336.2

参数与假数据中使用的参数有很大不同！我做错了什么？

【问题讨论】：

这可能不是问题的直接原因，但您似乎正在生成两个指数分布的加权平均值。但通常混合模型以w 的概率选择其中之一。
@KotaMori 谢谢。我不确定我是否在关注你。你的意思是我只需要估计lambda_1 和w？但是如何从中获得lambda_2？
不，我的意思是您生成假数据的代码可能不正确。请参阅下面的答案。

标签： r minimization mle

【解决方案1】：

原始代码容易出现警告和错误，因为参数可能很容易变为无效值。例如，我们需要w in [0, 1] 和lambda > 0。此外，如果a 大于一个数据点，则密度变为零，因此对数似然无限。

下面的代码使用了一些技巧来处理这些情况。

w 通过逻辑函数转换为范围 [0, 1]
lambda 通过指数函数转换为正值。
为处理零可能性案例的可能性增加了微小的价值。

此外，数据生成过程已更改，以便从具有给定概率w 的指数分布之一生成样本。

最后，由于n=500 的结果不稳定，所以增加了样本量。

biexpLL <- function(theta, y) {
  # define parameters
  w <- 1/(1+exp(-theta[1]))
  lambda_1 <- exp(theta[2])
  a <- theta[3]
  lambda_2 <- exp(theta[4])
  # likelihood function with dexp
  l <- w * dexp((y - a), rate = 1/lambda_1) + (1 - w) * dexp((y - a), rate = 1/lambda_2)
  - sum(log(l + 1e-9))
}
# Generate some fake data
w <- 0.7
n <- 5000
lambda_1 <- 2
lambda_2 <- 0.2
set.seed(45)
n1 <- round(n*w)
n2 <- n - n1
biexp_data <- c(rexp(n1, rate=1/lambda_1),
                rexp(n2, rate=1/lambda_2)) 
# Optimization
o <- optim(par=c(0.5,0.1,0.001,0.2),
           fn=biexpLL,
           y=biexp_data)

1/(1+exp(-o$par[1]))
exp(o$par[2])
o$par[3]
exp(o$par[4])

在我的环境中，我获得了以下内容。
结果似乎与模拟参数相当接近（请注意，交换了两个 lambda 值）。

> 1/(1+exp(-o$par[1]))
[1] 0.3458264
> exp(o$par[2])
[1] 0.1877655
> o$par[3]
[1] 3.738172e-05
> exp(o$par[4])
[1] 2.231844

请注意，对于这种混合模型，人们通常使用 EM 算法来优化似然性，而不是像这样直接优化。您可能也想看看它。

【讨论】：