[R-br] Dúvida em erro na função Filtro de Kalman

[Victor Eduardo]

Galera, estou com um problema que não estou entendendo o motivo do erro que
está na função filtro_kalman. Estou realizando um exercício que é para
construir uma função que realiza o filtro de kalman. Utilizei o modelo de
tendencia linear local no filtro de kalman:



#Função para gerar as matrizes do sistema para o
#Modelo de tendencia linear local

matriz_Tll<-function(fi,sig2_eta,sig2_eps,sig2_qui,n){

#Definindo as dimensões das matrizes do sistema

ar_Z<-array( ,dim=c(1,2,n))
ar_T<-array( ,dim=c(2,2,n))
ar_C<-array( ,dim=c(2,1,n))
ar_D<-array( ,dim=c(1,1,n))
ar_R<-array( ,dim=c(2,2,n))
ar_H<-array( ,dim=c(1,1,n))
ar_Q<-array( ,dim=c(2,2,n))

#Iniciando um 'for' para 'construir' as matrizes do sistema
for(t in 1:n){
ar_Z[ , ,t]<-matrix(c(1,0),nc=2)
ar_T[ , ,t]<-matrix(c(1,0,1,fi),nc=2)
ar_C[ , ,t]<-matrix(c(0,0))
ar_D[ , ,t]<-matrix(0)
ar_R[ , ,t]<-matrix(c(1,0,0,1),nc=2)
ar_H[ , ,t]<-matrix(sig2_eps,nc=1)
ar_Q[ , ,t]<-matrix(c(sig2_eta,0,0,sig2_qui),nc=2)
}
return(list(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q))
}

fi<-0.5 #com |fi|<1 temos um modelo de tendencia amortecida
sig2_eta<-4
sig2_eps<-6
sig2_qui<-10
n<-100 #número pequeno apenas para teste

mat_sist<-matriz_Tll(fi,sig2_eta,sig2_eps,sig2_qui,n)

ar_Z<-mat_sist[[1]]
ar_T<-mat_sist[[2]]
ar_C<-mat_sist[[3]]
ar_D<-mat_sist[[4]]
ar_R<-mat_sist[[5]]
ar_H<-mat_sist[[6]]
ar_Q<-mat_sist[[7]]


#Função para modelo em espaço de estado em geral

gerador_EE<-function(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1){
n<-dim(ar_Z)[3] #nº de matrizes do array Z ou D ou H
p<-dim(ar_H)[1] #nº linhas de H, D ou Z
m<-dim(ar_T)[1] #nº linhas de T, C ou R
r<-dim(ar_Q)[1] #nº linhas de Q
ar_Y<-array(,dim=c(p,1,n))
eps<-array(,dim=c(p,1,n))
eta<-array(,dim=c(r,1,n))
alfa<-array(,dim=c(m,1,n+1))
alfa[,,1]<-a1+chol(p1)%*%matrix(rnorm(m,0,1),nc=1)
for(t in 1:n){
eps[,,t]<-0+chol(ar_H[,,t])%*%matrix(rnorm(p,0,1),nc=1)
ar_Y[,,t]<-ar_Z[,,t]%*%alfa[,,t]+ar_D[,,t]+eps[,,t]
eta[,,t]<-0+chol(ar_Q[,,t])%*%matrix(rnorm(r,0,1),nc=1)
alfa[,,t+1]<-ar_T[,,t]%*%alfa[,,t]+ar_C[,,t]+ar_R[,,t]%*%eta[,,t]
}
return(ar_Y)
}

a1<-matrix(1:2,nc=1)
p1<-matrix(c(3,1,1,2),nc=2)
geracao<-gerador_EE(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1)

#Gráfico dos valores gerados para Y
plot(geracao,type='l')

#Filtro de Kalman

install.packages('MASS')
require(MASS)

Filtro_Kalman<-function(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1){

n<-dim(ar_Z)[3] #nº de matrizes do array Z ou D ou H
p<-dim(ar_H)[1] #nº linhas de H, D ou Z
m<-dim(ar_T)[1] #nº linhas de T, C ou R
r<-dim(ar_Q)[1] #nº linhas de Q

ar_Y<-array(,dim=c(p,1,n))
ar_U<-array(,dim=c(p,1,n))
ar_F<-array(,dim=c(p,p,n))
ar_AA<-array(,dim=c(m,1,n))
ar_PP<-array(,dim=c(m,m,n))
ar_A<-array(,dim=c(m,1,n+1))
ar_P<-array(,dim=c(m,m,n+1))

ar_A[,,1]<-a1
ar_P[,,1]<-p1
 ar_Y<-gerador_EE(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1)

for(t in 1:2){
ar_U[,,t]<-ar_Y[,,t]-ar_Z[,,t,drop=FALSE]%*%ar_A[,,t,drop=FALSE]-ar_D[,,t,drop=FALSE]
ar_F[,,t]<-ar_Z[,,t,drop=FALSE]%*%ar_P[,,t,drop=FALSE]%*%t(ar_Z[,,t,drop=FALSE])+ar_H[,,t,drop=FALSE]
ar_AA[,,t]<-ar_A[,,t,drop=FALSE]+ar_P[,,t,drop=FALSE]%*%ar_Z[,,t,drop=FALSE]%*%ginv(ar_F[,,t,drop=FALSE])%*%ar_U[,,t,drop=FALSE]
ar_PP[,,t]<-ar_P[,,t,drop=FALSE]-ar_P[,,t,drop=FALSE]%*%t(ar_Z[,,t,drop=FALSE])%*%ginv(ar_F[,,t,drop=FALSE])%*%ar_Z[,,t,drop=FALSE]%*%ar_P[,,t]
ar_A[,,t+1]<-ar_T[,,t,drop=FALSE]%*%ar_AA[,,t,drop=FALSE]+ar_C[,,t,drop=FALSE]
ar_P[,,t+1]<-ar_T[,,t,drop=FALSE]%*%ar_PP[,,t,drop=FALSE]%*%t(ar_PP[,,t,drop=FALSE])+ar_R[,,t,drop=FALSE]%*%ar_Q[,,t,drop=FALSE]%*%t(ar_R[,,t,drop=FALSE])
 }
return(list(ar_A,ar_P))
}

Filtro_Kalman(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1)
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120514/39d2ac20/attachment.html>