5

次元2x2の行列パラメーターの可能性を最大化しようとしています。尤度関数は、尤度も関数であるいくつかの固定行列パラメーターを渡す必要があります。Yで示されるデータと、共分散行列Sigma.star（下三角行列として通過します）が計算に必要ですが、これらを固定して、これに対して最適な関数を実行したいと思います。 Aを最適化しようとするコード

私の問題は、行列代数に使用しているオブジェクト内の何かを最適化しているという事実から、optimが誤っているように見えることです。すべての小さな計算をプログラミングせずにそれを機能させる方法はありますか？

具体的なエラーは次のとおりです。

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

しかし、クロネッカー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 に答える