R二元正态中的MLE答案

【问题标题】：MLE in R bivariate normalR二元正态中的MLE
【发布时间】：2018-04-22 19:31:18
【问题描述】：

我遇到了问题，可能是由于我的编码错误。我想通过算法对双变量正态样本执行 MLE：

require(MASS)
require(tmvtnorm)
require(BB)
require(matrixcalc)

mu = c(0,0)
covmat = matrix(c(9,-2,-2,4),2,2)

set.seed(357)

sample = list(rmvnorm(10,mu,covmat),
          rmvnorm(100,mu,covmat),
          rmvnorm(500,mu,covmat))

我设置了上面的示例。我将负对数似然定义如下；

neg_ll = function(mean_vec,cov_mat)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=mean_vec, sigma=cov_mat, log=TRUE))
  return(-log_ll)
}

当我运行以下 MLE 代码时，它会返回错误；

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean_vec=xbar_vec, cov_mat=scov_mat))

优化错误(start, f, method = method, hessian = TRUE, ...) : (list) 对象不能被强制输入'double'

NLM 函数有效（尽管仅适用于平均估计，不适用于协方差矩阵）。 NLM 返回：

nlm(f = neg_ll,xbar_vec, var(sample[[1]]),hessian = T)
$minimum
[1] 35.01874

$estimate
[1] -0.4036168  0.4703263

$gradient
[1] 1.463718e-06 3.886669e-06

$hessian
          [,1]      [,2]
[1,]  2.934318 -1.049366
[2,] -1.049366  7.769508

$code
[1] 1

$iterations
[1] 0

如何获得所有参数的估计值？我应该怎么做才能使用 MLE 功能？

编辑：neg_ll 函数出现类型错误 mean=mean 被 mean = mean_vec 替换。尽管如此，问题仍然存在，nlm 已估计仅平均向量的输出。

【问题讨论】：

能否在代码中包含您使用的库？
对不起，我忘记了，现在已经添加了。

标签： r statistics mle nlm log-likelihood

【解决方案1】：

如果您查看mle 使用的optimize 函数的注释，它说pars 应该是一维的。这就是您收到错误的原因。如果你像这样重写你的代码：

neg_ll = function(mean1, mean2, cov11, cov12, cov21, cov22)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=c(mean1, mean2), 
                   sigma=matrix(c(cov11, cov12, cov21, cov22), 2, 2), log=TRUE))
  return(-log_ll)
}

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
             cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], 
             cov21 = scov_mat[2, 1], cov22 = scov_mat[2, 2]))

您可以克服错误，但由于协方差矩阵中的非对角性，它会抛出一个新错误。因此，由于 dmvnorm 需要一个对角矩阵，因此您只需要上三角形或下三角形中的 1 个，在本例中为 1 个元素，cov12 或 cov21。

所以代码应该是这样的：

neg_ll = function(mean1, mean2, cov11, cov12, cov22)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=c(mean1, mean2), 
                   sigma=matrix(c(cov11, cov12, cov12, cov22), 2, 2), log=TRUE))
  return(-log_ll)
}

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
             cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], cov22 = scov_mat[2, 2]))

它给了我输出：

Call:
mle(minuslogl = neg_ll, start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
    cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], cov22 = scov_mat[2, 
        2]))

Coefficients:
     mean1      mean2      cov11      cov12      cov22 
-0.4036168  0.4703262  3.2228188  0.4352799  1.2171644

【讨论】：