<br />
kalmanl <- function(F,H,A,Q,R,Y,Z0,P0){ #return only the negative log likelihood value</p>
<p> N = ncol(Y)<br />
Y_0=list()<br />
V=list()<br />
Z_1=list()<br />
P_1=list()<br />
K=list()<br />
Z_0=list()<br />
P_0=list()<br />
loglike=list()<br />
ln=list()<br />
sln=NULL</p>
<p> Z_0[[1]]=Z0<br />
P_0[[1]]=P0<br />
I <-diag(1,nrow(Z0),nrow(Z0))</p>
<p> #filtering<br />
for (i in 1:N) {<br />
#forecasting Y<br />
Y_0[[i]]=t(H)%*%Z_0[[i]]<br />
V[[i]]=t(H)%*%P_0[[i]]%*%H+R</p>
<p> #updating the inference about Zt</p>
<p> Z_1[[i]]=Z_0[[i]]+ P_0[[i]] %*% H %*% solve(V[[i]]) %*% (Y[,i]-t(H)%*%Z_0[[i]])</p>
<p> P_1[[i]]=P_0[[i]]- P_0[[i]] %*% H %*% solve(V[[i]]) %*% t(H)%*% P_0[[i]]</p>
<p> # Kalman Gain<br />
K[[i]]=F%*%P_0[[i]]%*%H%*%solve(V[[i]])</p>
<p> #producing a forecast of Z<br />
Z_0[[i+1]]=F%*%Z_1[[i]]+A # A is the const in the model</p>
<p> P_0[[i+1]]=F%*%P_1[[i]]%*%t(F)+Q</p>
<p> # return the loglikelihood </p>
<p> loglike[[i]]=-log(det(V[[i]]))-t(Y[,i]-Y_0[[i]])%*%solve(V[[i]])%*%(Y[,i]-Y_0[[i]]) </p>
<p> }</p>
<p> # log likelihood function<br />
ln=do.call(cbind,loglike)<br />
sln=-sum(ln)<br />
result=list(sln=sln) </p>
<p> return(result)<br />
}</p>
<p>
我用optim()来优化上面的likelihood function.就是 <bblatex> loglike=-log(det(V))-t(Y-Y_0])%*%solve(V)%*%(Y-Y_0) </bblatex>
但是优化到一半就停住了,因为矩阵V[
]是singular matrix,没有det和逆矩阵. 请问如何修改上面这个function可以让optim跳过错误的情况继续优化?
</p>
