如何解决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算法来优化可能性,而不是像这样直接优化。您可能还想看看它。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。