0

我想最小化函数 FlogV（使用多正态分布，Z 是数据矩阵 NxC；SIGMA 它是数据的 var-covariance 的方阵 CxC，R 是长度为 C 的向量）

FLogV <- function(P){

(here I define parameters, P, within R and SIGMA)

logC <- (C/2)*N*log(2*pi)+(1/2)*N*log(det(SIGMA))

SOMA.t <- 0

for (j in 1:N){
SOMA.t <- SOMA.t+sum(t(Z[j,]-R)%*%solve(SIGMA)%*%(Z[j,]-R))
}

MlogV <- logC + (1/2)*SOMA.t

return(MlogV)

}

minLogV <- optim(P,FLogV)

所有这些都是扩展代码的一部分，该代码已经过测试并且运行良好，除了最重要的事情：我无法优化，因为我收到此错误：“solve.default(SIGMA) 中的错误：系统在计算上是奇异的：倒数条件编号 = 3.57726e-55”</p>

如果我使用 ginv() 或 pseudoinverse() 或 qr.solve() 我会得到：“svd(X) 中的错误：'x' 中的无限或缺失值”</p>

问题是：如果我在错误消息之后取 SIGMA 矩阵，我可以求解（SIGMA），特征值都是正的，行列式非常小，但是正的 det(SIGMA) [1] 3.384674e-76

eigen(SIGMA)$values
[1] 0.066490265 0.024034173 0.018738777 0.015718562 0.013568884 0.013086845
….
[31] 0.002414433 0.002061556 0.001795105 0.001607811

我已经阅读了几篇关于像 SIGMA 之类的变化矩阵（接近单数）的论文，对数据规模和形式进行了几次转换，但我意识到，对于像示例这样的 34x34 矩阵，在 det(SIGMA) 接近 e-40 之后， R假设它像0并且计算失败；我也不能减少矩阵维数，也不能在我的函数校正算法中输入奇异矩阵，因为 R 不能使用像 optim 这样的优化函数来评估它。我非常感谢对这个问题的任何建议。在此先感谢，玛丽亚 D。

4

1 回答 1

1

从您的帖子中不清楚故障是来自det()还是solve()

如果它只是二次项中的解，您可能想尝试 的两个参数版本solve，它可以更稳定一些。solve(X,Y)是相同的solve(X) %*% Y

如果您可以使用 分解 sigma chol()，您将得到一个三角矩阵，使得 LL'=Sigma。行列式是对角线的乘积，你可以试试这个二次项：

crossprod( backsolve(L, Z[j,]-R))
于 2013-11-20T01:34:33.073 回答