如何在R中找到特定多元正态对数似然函数的最大似然值？

Question

我在R中优化多元正态对数似然时遇到了麻烦，如果有人有好的解决方案，请告诉我。具体来说，我似乎无法保持协方差矩阵正定并使参数处于合理范围内。

让我更全面地介绍一下这个问题。我基本上正在尝试使用MLE同时解决这两个回归方程：

$$ y_1 = \beta_1 + \beta_2 x_1 + \beta_3 x_2 \\ y_2 = \beta_4 + \beta_3 x_1 + \beta_5 x_2 $$

事实上，$\beta_3$同时出现在两个方程式中并不是一个错误。我通过最大化多元正态分布的似然性来解决这个问题，其中$Y = (y_1, y_2)^\top$的均值像上述回归方程那样进行参数化。

我已经附上了对数似然函数，因为我认为它应该是这样的。我通过重新创建必须为正的特征值和Cholesky分解来将方差协方差矩阵限制为正定。

mvrestricted_ll <- function(par, Y, X) {

  # Indices
  n <- nrow(X)
  nbetas <- (2 + 3 * (ncol(Y) - 1))

  # Extract parameters
  beta <- par[1:nbetas]
  eigvals <- exp(par[(nbetas + 1):(nbetas + ncol(Y))]) # constrain to be positive
  chole <- par[(nbetas + ncol(Y) + 1):(nbetas + ncol(Y) + ncol(Y)*(ncol(Y)+1)/2)]

  # Build Sigma from positive eigenvalues and cholesky (should be pos def)
  L <- diag(ncol(Y))
  L[lower.tri(L, diag=T)] <- chole
  Sigma <- diag(eigvals) + tcrossprod(L)

  # Linear predictor
  # Hard coded for 2x2 example for now
  mu <- cbind(beta[1] + beta[2]*X[,1] + beta[3]*X[,2],
              beta[4] + beta[3]*X[,1] + beta[5]*X[,2])

  yminmu <- Y - mu

  nlogs <- n * log(det(Sigma))

  invSigma <- solve(Sigma)

  meat <- yminmu %*% tcrossprod(invSigma, yminmu)

  return(- nlogs - sum(diag(meat)))
}

# Create fake data
n <- 1000
p <- 2
set.seed(20160201)
X <- matrix(rnorm(n*p), nrow = n)
set.seed(20160201)
Y <- matrix(rnorm(n*p), nrow = n)

# Initialize parameters
initpars <- c(rep(0, (2 + 3 * (ncol(Y) - 1)) + ncol(Y) + ncol(Y)*(ncol(Y)+1)/2))
# Optimize fails with BFGS
optim(par = initpars, fn = mvrestricted_ll, X=X, Y=Y, method = "BFGS")
# Optim does not converge with Nelder-mead, if you up the maxits it also fails
optim(par = initpars, fn = mvrestricted_ll, X=X, Y=Y)

非常感谢您的帮助。

编辑：需要注意的是，仅在参数中让Sigma成为向量，然后在它不是正定时返回一个非常大的值也不起作用。

- stich09

Answer 1

我不确定代码/答案是否正确，但是

invSigma <- try(solve(Sigma))
if (inherits(invSigma, "try-error")) return(NA)

并运行

optim(par = initpars, fn = mvrestricted_ll, X=X, Y=Y,
      control = list(maxit = 1e5))

我已经取得了一些进展，将收敛代码提高到10（退化的Nelder-Mead单纯形）。

$par
 [1]  1.361612e+01  4.674349e+01 -3.050170e+01  3.305013e+01  6.731194e+01
 [6] -3.117192e+01 -5.408598e+00 -6.326897e-07 -1.987449e+01 -1.795924e+01

$value
[1] -1.529013e+19

$counts
function gradient 
    1219       NA 

$convergence
[1] 10

我猜想一个真正的解决方案可能需要更仔细地查看代码，以确定它是否真的在做你认为它在做的事情（抱歉）；了解为什么会出现 solve() 错误可能是一个很好的第一步。您可以通过将 cat(par, "\n") 放在函数的第一行，然后运行 没有 try/NA-return 代码来解决此问题。这样可以使您隔离出导致错误的示例数据集，然后您可以逐行（使用 debug() 或手动方式）逐行分析您的代码以了解发生了什么。

- Ben Bolker

Answer 2

您可以考虑使用以下方法：

library(DEoptim)

fn <- function(par, mat_X, mat_Y)
{
  X <- mat_X
  Y <- mat_Y
  n <- nrow(X)
  nbetas <- (2 + 3 * (ncol(Y) - 1))
  beta <- par[1 : nbetas]
  eigvals <- exp(par[(nbetas + 1) : (nbetas + ncol(Y))]) 
  chole <- par[(nbetas + ncol(Y) + 1) : (nbetas + ncol(Y) + ncol(Y) * (ncol(Y) + 1) / 2)]
  L <- diag(ncol(Y))
  L[lower.tri(L, diag = TRUE)] <- chole
  
  Sigma <- tryCatch(diag(eigvals) + tcrossprod(L), error = function(e) NA)
  
  if(is.null(dim(Sigma)))
  {
    return(10 ^ 30)
    
  }else
  {
    mu <- cbind(beta[1] + beta[2] * X[,1] + beta[3] * X[,2],
                beta[4] + beta[3] * X[,1] + beta[5] * X[,2])
    
    yminmu <- Y - mu
    nlogs <- n * log(det(Sigma))
    invSigma <- tryCatch(solve(Sigma), error = function(e) NA)
    
    if(is.null(dim(invSigma)))
    {
      return(10 ^ 30)
      
    }else
    {
      meat <- yminmu %*% tcrossprod(invSigma, yminmu)
      log_Lik <- - nlogs - sum(diag(meat))
      
      if(is.na(log_Lik) | is.nan(log_Lik) | is.infinite(log_Lik))
      {
        return(10 ^ 30)
        
      }else
      {
        return(-log_Lik)
      }
    }
  }
}

n <- 1000
p <- 2

set.seed(20160201)
mat_X <- matrix(rnorm(n * p), nrow = n)

set.seed(2436537)
mat_Y <- matrix(rnorm(n * p), nrow = n)

lower <- rep(-10, 10)
upper <- rep(10, 10)
DEoptim(fn = fn, lower = lower, upper = upper, 
        control = list(itermax = 10000, parallelType = 1), mat_X = mat_X, mat_Y = mat_Y)

- Emmanuel Hamel