Changes made in revision 6 of KalmanEM.r
--- /tmp/cvx_fYRpTv 2012-02-07 19:48:37.000000000 -0800
+++ /tmp/cvx_xUwD3l 2012-02-07 19:48:37.000000000 -0800
@@ -1,7 +1,8 @@
KalmanEM = function(y, whichPop=NA, Binit=NA, Uinit=0, Qinit=0.05, Ainit=0, Rinit=0.05, x00=NA, V00init=10, Zinit=NA, max.iter=5000, varcov.Q="unconstrained", varcov.R="diagonal", B.constraint="identity", U.constraint="unconstrained", Z.constraint="missingvalues", A.groups=NA, U.groups=NA, Q.groups=NA, R.groups=NA, x00.groups=NA, B.groups=NA, U.bounds=c(-1,1), logQ.bounds = c(log(1.0e-05),log(1.0)), logR.bounds = c(log(1.0e-05),log(1.0)), x00prior=FALSE, estInteractions=FALSE, MonteCarloInit = FALSE, numInits = 500, numInitSteps = 10, tol=0.01, silent=FALSE, debugit=FALSE)
{
#Version
-vers = "2.68" #EEH Sept 24, 2009
+vers = "2.69" #EEH Oct 4, 2009
+#KalmanEM written by Eli Holmes and Eric Ward
#comment on x00prior=FALSE; x00prior=TRUE case is not programmed in yet. In that case x00 fixed, not estimated.
#x00prior=FALSE means that we treat x00 as an unknown fixed parameter that is estimated. So P_0=V00=0, but we don't want to set
#V00=0 because then the EM algorithm won't converge. We want V00 big, and then at the end, return the likelihood with V00=0
@@ -510,11 +511,11 @@
cat(paste("Finished in ",iter," interations. Max.iter was ",max.iter,".\n",sep=""))
}
# confidence intervals based on state std errors, see caption of Fig 6.3 (p337) Shumway & Stoffer
-if(m==1) se = kf$VtT[,,1:numYrs]
-if(m > 1) {
- se = matrix(0,nrow=m,ncol=numYrs)
- for(i in 1:numYrs) se[,i] = t(sqrt(takediag(kf$VtT[,,i])))
-}
+#if(m==1) states.se = sqrt(kf$VtT[,,1:numYrs])
+#if(m > 1) {
+ states.se = matrix(0,nrow=m,ncol=numYrs)
+ for(i in 1:numYrs) states.se[,i] = t(sqrt(takediag(kf$VtT[,,i])))
+#}
# This calculates the correct likelihood when x00 is estimated, that is treated as fixed but unknown
# Since it is fixed, V00=0 but we can't run the EM algorithm with V00=0 or it would not converge
if(x00prior==FALSE) {
@@ -563,7 +564,7 @@
colnames(B) = state.names
rownames(kf$xtT) = state.names
-return(list(states=kf$xtT,A=A,B=B,Q=Q,R=R,U=U,Z=Z,x00=x00,V00=V00,V0T=kf$V0T,Kt=kf$Kt,Innov=kf$Innov,Sigma=kf$Sigma,m=m, n=n,
+return(list(states=kf$xtT,states.se=states.se,A=A,B=B,Q=Q,R=R,U=U,Z=Z,x00=x00,V00=V00,V0T=kf$V0T,Kt=kf$Kt,Innov=kf$Innov,Sigma=kf$Sigma,m=m, n=n,
numYrs=numYrs,iter=iter,loglike=loglike, AIC=AICval, AICc=AICc, iter.record=iter.record,
whichPop=whichPop, R.groups=R.groups, Q.groups=Q.groups, U.groups=U.groups, varcov.R=varcov.R, varcov.Q=varcov.Q,
estInteractions=estInteractions, x00prior=x00prior, data=orig.data, kf.data=y, include=include, version=vers))