r - 如何在R中最大化函数内部使用矩阵乘法来获得最佳效果

Question

我正在尝试使尺寸为2x2的Matrix参数的可能性最大化。似然函数需要传入两个固定矩阵参数，似然也是函数。表示Y的数据和协方差矩阵Sigma.star(我作为下三角矩阵传递)是计算所必需的，但我想保持这些固定不变，并在此基础上运行优化函数尝试优化A的代码

我的问题是，它似乎正在优化我正在用于矩阵代数的对象内部的某些东西，这似乎是一种错误。有什么方法可以使它工作而无需编写每一个小小的计算？

具体错误是:

Error in diag(1, nrow = (m^2)) - A %x% A : non-conformable arrays

但是A kronecker A应该是一个m ^ 2 x m ^ 2矩阵，就像恒等式一样。

代码:

library(MCMCpack)
library(mvtnorm)
set.seed(1000)


Likelihood.orig<-function(A, Y, Sigma.star){
Sigma<-xpnd(Sigma.star)
n<-nrow(Y)
if(is.vector(A)==TRUE){
  A<-as.matrix(A, nrow=nrow(Sigma), ncol=ncol(Sigma))
}
m<-nrow(A)
V<-matrix(solve(diag(1, nrow=(m^2))-A%x%A)%*%as.vector(Sigma), nrow=m, ncol=m)
temp1<- (-.5)*log(abs(det(V)))
temp2<- (-(n-1)/2)*log(abs(det(Sigma)))
temp3<- t(Y[,1, drop=FALSE]) %*% (solve(V)) %*% Y[,1, drop=FALSE]
terms<- numeric(n-1)
for(i in 2:n){
    terms[i-1]<- t(Y[,i, drop=FALSE] - A %*%Y[,i-1, drop=FALSE]) %*% (solve(Sigma)) %*% (Y[,i] - A %*%Y[,i-1])
}
return(temp1+temp2-.5*(temp3+sum(terms)))   
}




Generate.Y<-function(n, A, Sigma){
m<-nrow(A)
Y<-matrix(0, nrow=m, ncol=n)
V<-matrix(solve(diag(1, nrow=m^2)-A%x%A)%*%as.vector(Sigma), nrow=m, ncol=m)
Y[,1]<-rmvnorm(1, numeric(nrow(A)), V)
for(i in 2:n){
    Y[,i]<-A%*%Y[,i-1, drop=FALSE]+t(rmvnorm(1, mean = numeric(m), sigma = Sigma))
    }
return(Y)
}


n<-500
A.true<-matrix(c(.8, .3, 0, .5), nrow=2, ncol=2)
Sigma<-matrix(c(1, 0, 0, .5), nrow=2, ncol=2)
Y<-matrix(0, nrow=2, ncol=n)
Y<-Generate.Y(n, A.true, Sigma)
m=nrow(Y)
lower.Sigma<-vech(Sigma)



optim(par=c(1, 0, 0, 1), fn=Likelihood.orig, method="Nelder-Mead", 
  control=list(maxit=500, fnscale=-1), Sigma.star=lower.Sigma, Y=Y)

Answer 1

您的方法是正确的，即让optim在向量上进行优化，并且仅将该向量转换为要最大化的函数内的矩阵。

但是，您需要使用matrix而不是as.matrix来创建该矩阵。查看以下内容之间的区别:

as.matrix(1:4, nrow=2, ncol=2)  # wrong tool
#      [,1]
# [1,]    1
# [2,]    2
# [3,]    3
# [4,]    4

和

matrix(1:4, nrow=2, ncol=2)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

对于此类问题，我强烈建议您学习R调试工具( browser， debug， debugonce等)。有关示例，请参见 General suggestions for debugging in R。