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))

Sculpin 0.2 | xhtml | problems or comments? | report bugs