R双变量正常中的MLE

Question

我遇到了问题，可能是因为我的编码错误。我想通过算法对二元正态样本执行 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 代码时 returns 一个错误；

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))

Error in optim(start, f, method = method, hessian = TRUE, ...) : (list) object cannot be coerced to type 'double'

NLM 函数有效（尽管仅适用于均值估计，不适用于协方差矩阵）。 NLM returns:

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 估计平均向量的输出仅。

Answer 1

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

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

R双变量正常中的MLE

MLE in R bivariate normal

statistics

r

nlm

mle

log-likelihood