1

我遇到了问题，可能是由于我的编码错误。我想通过算法对双变量正态样本执行 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))

optim(start, f, method = method, hessian = TRUE, ...) 中的错误：（列表）对象不能被强制输入“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仅估计了平均向量的输出。

4

1 回答 1

1

如果你看一下optimize函数的注释mle，它说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
于 2017-11-09T19:28:04.187 回答