5

我试图最大化维度 2x2 的矩阵参数的可能性。似然函数需要传入几个固定矩阵参数，似然函数也是它们的函数。表示为 Y 的数据和协方差矩阵 Sigma.star（我将其作为下三角矩阵传递）是计算所必需的，但我想保持这些数据固定并在我的试图优化 A 的代码

我的问题是 optim 似乎是错误的，因为它正在优化我用于矩阵代数的对象内部的某些内容。有没有什么方法可以让它工作而不需要对每一个小计算进行编程？

具体错误是：

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

但是 A kronecker A 应该是一个 m^2 xm^2 矩阵，就像恒等式一样......</p>

代码：

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

1 回答 1

7

您的方法是正确的，即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 中调试的一般建议。

于 2012-11-15T00:02:14.897 回答