Comment optimiser le fonctionnement de la multiplication matricielle dans la fonction à maximiser dans R

Question

Je tente de maximiser la probabilité d'un paramètre Matrix de dimension 2x2. La fonction de vraisemblance doit transmettre un couple de paramètres matriciels fixes dont la vraisemblance est également une fonction. Les données, notées Y, et une matrice de covariance, Sigma.star (que je traverse en tant que matrice triangulaire inférieure), sont nécessaires pour le calcul, mais je voudrais les garder fixes et exécuter une fonction optim sur ceci, dans mon code essayant d'optimiser A

Mon problème est que optim semble se tromper du fait qu'il optimise quelque chose à l'intérieur d'un objet que j'utilise pour l'algèbre matricielle. Y a-t-il un moyen de le faire fonctionner sans programmer chaque petit calcul?

L'erreur spécifique est:

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

Mais un kronecker A devrait être un m^2 xm^2 matrice comme l'identité ...

code:

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

Votre approche est correct, c'est-à-dire, faites optim optimiser sur un vecteur, et ne faites que transformer ce vecteur en une matrice dans la fonction que vous essayez de maximiser. Toutefois, vous devez utiliser matrix et non as.matrix pour créer cette matrice. Voyez la différence entre:

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

et

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

Pour les problèmes de ce type, je vous recommande vivement d'apprendre les outils de débogage R (browser, debug, debugonce, etc.). Voir General suggestions for debugging in R pour des exemples.