Changes made in revision 3 of KalmanEM.r
--- /tmp/cvx_7Klr4C 2012-02-07 19:52:25.000000000 -0800
+++ /tmp/cvx_NydRjr 2012-02-07 19:52:25.000000000 -0800
@@ -1,5 +1,8 @@
-KalmanEM = function(y, whichPop=NA, Binit=NA, Uinit=0, Qinit=0.05, Ainit=0, Rinit=0.05, x00=NA, V00init=10, max.iter=5000, varcov.Q="unconstrained", varcov.R="diagonal", B.constraint="identity", U.constraint="unconstrained", 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)
-{ #comment on x00prior=FALSE; x00prior=TRUE case is not programmed in yet. In that case x00 fixed, not estimated.
+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.67" #EEH July 17, 2009
+ #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
# Code written by Eric Ward (EW) and Eli Holmes (EH)
@@ -8,8 +11,6 @@
# The model needs to specified in a log-linear form; for population count data, this means log-transforming the data
# Will deal with missing values (marked as -99), but this slows computation
-
-
############## Things that the function needs passed in
###########################################################
# y is the LOG TRANSFORMED data with -99 for MISSING VALS; each column is a separate observation time series; time is down the rox
@@ -82,7 +83,7 @@
############ Error check that the model is legal or consistent
###########################################################
-if(FailedErrorCheck(y, whichPop, U.groups, x00.groups, A.groups, R.groups, Q.groups, B.groups, varcov.Q, varcov.R, B.constraint, Uinit, Qinit, Ainit, Rinit, Binit, x00, V00init,estInteractions))
+if(FailedErrorCheck(y, whichPop, U.groups, x00.groups, A.groups, R.groups, Q.groups, B.groups, varcov.Q, varcov.R, B.constraint, Uinit, Qinit, Ainit, Rinit, Binit, x00, V00init,estInteractions, silent))
stop("Model is illegal or not interally consistent")
############# Set the initial conditions and make them the correct dimensions
@@ -114,8 +115,11 @@
############# Construct the various design matrices
# Z relates observation time series to state processes
-Z = matrix(0,n,m) # this will be filled 0s and 1s...if each site is unique, Z = diag(n)
-for(i in 1:m) Z[which(whichPop==i),i] <- 1
+if(Z.constraint=="fixed") Z = Zinit #the user passes this in
+else { #Z is constructed and Zinit is not used
+ Z = matrix(0,n,m) # this will be filled 0s and 1s...if each site is unique, Z = diag(n)
+ for(i in 1:m) Z[which(whichPop==i),i] <- 1
+ }
# Set up the grouping structure among different state processes and the different observation time se
A.numGroups <- max(A.groups)
@@ -179,16 +183,19 @@
cat(modelStruc)
}
-#############Set up some matrices used for handling missing values
+#############Adapt the Z matrix for handling missing values
# include is a n x m matrix to handle missing values; 1 matrix for each year;
# if a data point in time series j in yr i is missing, element j in the diagonal of the include matrix for year i is set to 0
###########################################################
-include <- array(0, dim=c(n,m,numYrs))
-for(i in 1:numYrs)
- include[,,i] <- makediag(ifelse(y[,i]!=-99,1,0),nrow=n)%*%Z #include = 1 when data exists, 0 otherwise
-miss <- ifelse(y!=-99,0,1) # 0=observed, 1=missing
-notmiss <- ifelse(y!=-99,1,0) # 1=observed, 0=missing
-y[which(y == -99)] <- 0 # Clear out -99s in data; include has been set up so that these values will be ignored
+if(Z.constraint=="fixed") include=Zinit
+else {
+ include <- array(0, dim=c(n,m,numYrs))
+ for(i in 1:numYrs)
+ include[,,i] <- makediag(ifelse(y[,i]!=-99,1,0),nrow=n)%*%Z #include = 1 when data exists, 0 otherwise
+ miss <- ifelse(y!=-99,0,1) # 0=observed, 1=missing
+ notmiss <- ifelse(y!=-99,1,0) # 1=observed, 0=missing
+ y[which(y == -99)] <- 0 # Clear out -99s in data; include has been set up so that these values will be ignored
+}
#############Start the EM algorithm which is going to keep updating until
# it converges
@@ -230,10 +237,12 @@
# this determines # of iterations, if in MCInit, then iter = numInitSteps otherwise max.iter
numIter = ifelse(loop <= numInits,numInitSteps,max.iter)
if(MonteCarloInit == TRUE) { # draw random initial conditions / parameter values
- U = array(runif(m,U.bounds[1],U.bounds[2]),dim=c(m,1))
+ if(U.constraint!="fixed") U = array(runif(m,U.bounds[1],U.bounds[2]),dim=c(m,1))
Q = makediag(exp(runif(m,logQ.bounds[1],logQ.bounds[2])),nrow=m)
R = makediag(exp(runif(n,logR.bounds[1],logR.bounds[2])),nrow=n)
- x00 = array(runif(m,0.75*x00init,1.25*x00init),dim=c(m,1))
+ x.lo = ifelse(x00init>0,0.75*x00init,1.25*x00init)
+ x.hi = ifelse(x00init>0,1.25*x00init,0.75*x00init)
+ x00 = array(runif(m,x.lo,x.hi),dim=c(m,1))
if(drawProgressBar==TRUE) prev = progressBar(loop/numInits,prev)
}
# Some initial start draws of the params are so awful that loglike will go down. We want to discard these cases.
@@ -296,6 +305,19 @@
} # end for subpop loop
A <- as.vector(c(ZA%*%(colSums(makediag(A,nrow=n)%*%ZA) * 1/colSums(ZA))) )
} # end if
+
+ ################
+ # EH: If we treat x00 as fixed and unknown, we estimate it.
+ # Get new x00, V00 cannot be estimated so it is fixed
+ # S&S Eqn 4.78
+ # If x00 is a known prior, we don't estimate.
+ ################################################################
+ if(x00prior==FALSE) {
+ x00 = c(kf$x0T)
+ x00 = as.vector(c(Zx00%*%(colSums(makediag(x00,nrow=m)%*%Zx00) * 1/colSums(Zx00))) ) #impose grouping
+ x00 = array(x00, dim=c(m,1))
+ }
+
################
# Get new U subject to its constraints (update of Q and B will use this)
# See notes by EH and Rich Hinrichsen for derivation
@@ -364,7 +386,16 @@
# S10 <- kf$Vtt1T[,,1] + (kf$xtT[,1]-U)%*%t(kf$x0T);
# I switch on 7/22/08 to this since it seems less sensitive to V00 and finds Q with max L with lower tol setting
# with S&S it climbs Q at more slowly and cvg hits tol before max is reached
+ # put x00 update earlier
+
+ # S00 <- V00 + x00%*%t(x00)
+ # S11 <- kf$VtT[,,1] + (kf$xtT[,1]-U)%*%t(kf$xtT[,1]-U)
+ # S10 <- kf$Vtt1T[,,1] + (kf$xtT[,1]-U)%*%t(x00);
+
################################################################
+ iter.model=list( U=U, Q=Q, R=R, A=A, B=B, x00=x00, V00=V00, include=include, kf.data=y )
+ kf <- Kfilter(iter.model) #call kalman filter/smoother to get x(t)|y(1:T) etc estimates
+ loglike.new = kf$loglik
S00 = 0; S11 = 0; S10 = 0; X1 = 0; X0 = 0
for (i in 2:numYrs) {
S00 = S00 + (kf$VtT[,,i-1] + kf$xtT[,i-1]%*%t(kf$xtT[,i-1])); #sum 2:T E(xt1T%*%t(xt1T))
@@ -375,16 +406,6 @@
}
################
- # EW: Get new B subject to its constraints
- ################################################################
- if(B.constraint!="fixed" & B.constraint!="identity")
- B_unconstrained = (S10-U%*%t(X0))%*%chol2inv(chol(S00))
- if(B.constraint == "unconstrained") B <- B_unconstrained
- if(B.constraint == "diagonal") { #allows shared B's along diagonal
- B=makediag(takediag(S10-U%*%t(X0))/takediag(S00))
- B <- makediag(c(ZB%*%(colSums(B%*%ZB) * 1/colSums(ZB))),nrow=m)
- }
- ################
# Get new Q subject to its constraints
################################################################
Q_unconstrained = (S11 - B%*%t(S10) - S10%*%t(B) + B%*%S00%*%t(B)
@@ -394,9 +415,8 @@
#The inv functions don't work well when the diagonal elements get lower than the machine tol, so put a min on those
if(length(Q_unconstrained)==1)
Q_unconstrained <- max(sqrt(.Machine$double.eps),Q_unconstrained)
- else
- diag(Q_unconstrained)[which(diag(Q_unconstrained)