R：使用 optim 的指数混合的最大似然估计

Question

我正在尝试使用 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

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

Answer 1

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

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

w被logistic函数

[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 算法来优化似然，而不是像这样直接优化。你可能也想看看它。

R：使用 optim 的指数混合的最大似然估计

R: Maximum Likelihood Estimation of a exponential mixture using optim

r

minimization

mle