如何在函数内部使用矩阵乘法进行优化,以便在R中最大化

Question

我试图最大化尺寸为2x2的Matrix参数的可能性.似然函数需要传递几个固定矩阵参数,其可能性也是函数.表示为Y的数据和协方差矩阵Sigma.star(我将其作为下三角矩阵传递)是计算所必需的,但我想保留这些数据并对此运行优化函数,在我的代码试图优化A.

我的问题是,似乎错误是因为它正在优化我用于矩阵代数的对象内部的事实.是否有一些方法可以使它工作而无需编程每一个小计算？

具体错误是:

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

但是一个kronecker A应该是一个m ^ 2 xm ^ 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等).有关示例,请参阅R中的调试常规建议.