在 r 中编码似然和对数似然函数以执行优化

Question

我正在写一篇论文，要求我找到 Gumbel 类型 I 的 MLE 双变量指数分布。我已经证明了似然函数和对数似然函数 likelihood and log-likelihood 但我正在努力在 r 中实现它以使用 Optim 函数执行优化。我的代码生成 NA 值。以下是我的代码。

# likelihood function of x
likelihood.x = function(params, data) {
  lambda1 = params[1]
  lambda2 = params[2]
  theta = params[3]
  A = (1 - theta) * (lambda1 * lambda2)
  B = theta * (lambda1 ^ 2) * lambda2 * data$X1
  C = theta * lambda1 * (lambda2 ^ 2) * data$X2
  D = (theta ^ 2) * (lambda1 ^ 2) * (lambda2 ^ 2) * data$X1 * data$X2
  E = (lambda1 * data$X1) + (lambda2 * data$X2) + (theta * lambda1 * lambda2 * data$X1 * data$X2)
  f = sum(log(A + B + C + D)) - sum(E)
  return(exp(f))
}

# Log-likelihood function of x
log.likelihood.x = function(params, data){
  lambda1 = params[1]
  lambda2 = params[2]
  theta = params[3]
  A = (1 - theta) * (lambda1 * lambda2)
  B = theta * (lambda1 ^ 2) * lambda2 * data$X1
  C = theta * lambda1 * (lambda2 ^ 2) * data$X2
  D = (theta ^ 2) * (lambda1 ^ 2) * (lambda2 ^ 2) * data$X1 * data$X2
  E = (lambda1 * data$X1) + (lambda2 * data$X2) + (theta * lambda1 * lambda2 * data$X1 * data$X2)
  f = sum(log(A + B + C + D)) - sum(E)
  return(-f)
}

这是生成数据的函数

# Simulating data
rGBVE = function(n, lambda1, lambda2, theta) {
  x1 = rexp(n, lambda1)
  lambda12 = lambda1 * lambda2
  pprod = lambda12 * theta
  C = exp(lambda1 * x1)
  A = (lambda12 - pprod + pprod * lambda1 * x1) / C
  B = (pprod * lambda2 + pprod ^ 2 * x1) / C
  D = lambda2 + pprod * x1
  wExp = A / D
  wGamma = B / D ^ 2
  data.frame(x1, x2 = rgamma(n, (runif(n) > wExp / (wExp + wGamma)) + 1, D))
}

data = rGBVE(n=100, lambda1 = 1.2, lambda2 = 1.4, theta = 0.5)
colnames(data) = c("X1", "X2")

我的目标是使用 r 中的 Optim() 找到 lambda1、lambda2 和 theta 的 MLE。

请帮助我在 r 中实现似然函数和对数似然函数。谢谢。

Answer 1

您的问题似乎与警告消息有关

In log(A+B+C+D): NaNs produced

此类警告通常无害 — 它只是意味着优化算法在某个地方尝试了一组参数，但违反了条件 A+B+C+D ≥ 0。由于这些是相当复杂的表达式，因此需要花费一些精力来弄清楚如何约束参数（或重新参数化函数，例如在对数刻度上拟合某些参数）以避免警告，但猜测保留参数 non-negative 会有所帮助，我们可以尝试使用 L-BFGS-B 算法（这是 optim() 中唯一可用的允许多维有界优化的算法）。

r1 <- optim(par = c(1,2,1),
      fn = log.likelihood.x,
      dat  = data)

r2 <- optim(par = c(1,2,1),
      fn = log.likelihood.x,
      lower = rep(0,3),
      method = "L-BFGS-B",
      dat  = data)

第二个不生成警告，结果接近（如果不相同）：

all.equal(r1$par, r2$par)
## "Mean relative difference: 0.0001451953"

您可能想要使用 bbmle，它具有一些用于可能性建模的附加功能：

library(bbmle)
fwrap <- function(x) log.likelihood.x(x, dat = data)
parnames(fwrap) <- c("lambda1", "lambda2", "theta")
m1 <- mle2(fwrap, start = c(lambda1 = 1, lambda2 = 2, theta = 1), vecpar = TRUE,
           method = "L-BFGS-B", lower = c(0, 0, -0.5))
pp <- profile(m1)
plot(pp)
confint(pp)
confint(m1, method = "quad")

在 r 中编码似然和对数似然函数以执行优化

Coding likelihood and log-likelihood function in r to perform optimization

optimization

r

exponential

mle