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

问题描述

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

这是代码

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 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算法来优化可能性，而不是像这样直接优化。您可能还想看看它。