KalmanEM = function(y, whichPop=NA, Uinit=NA, Qinit=NA, Ainit=NA, Rinit=NA, x00=NA, V00init=NA, max.iter=5000, varcov.Q="unconstrained", varcov.R="diagonal", U.groups=NA, Q.groups=NA, R.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)), estInteractions=FALSE, MonteCarloInit = FALSE, numInits = 500, numInitSteps = 10, tol=0.01, silent=FALSE)
{
# Code written by Eric Ward (EW) and Eli Holmes (EH)
# adapted from Matlab code by Eli Holmes/Rich Hinrichsen and R code from Shumway and Stoffer (2006)
# A description of the algorithm (for the univariate case) is in "An EM algorithm.pdf"
# July 4, 08, EH added estimation A bias term
# July 7, 08 EW added arguments whichPop (nx1), R.groups (nx1), U.groups (mx1) and Q.groups (mx1)
# July 12, 08 EH cleaned up error checking
# July 19, 08 EW added a few catches to trap cases where A is never estimated
# Aug 06, 08 EW added names, so that strings may be passed in for groups / sites / species (names also display on output)
# Nov 08 EW added bootstrapping functions
############## 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
y <- t(y) #for consistency, all data files are assumed to have time going down rows; this algorithm wants time across cols
# Initial conditions for the algorithm are passed in as vectors Uinit,Qinit,Rinit,x00
# max.iter is the number of iterations to do
# tol is the tolerence for checking if convergence has happened; lower if convergence happens in just a few iterations
# whichPop, a 1xn vector assigning observation time series to state processes
# varcov.Q and varcov.R specify the structures allowed for Q & R: "unconstrained", "diagonal" or "equalvarcov"
# R.groups is a nx1 vector specifying which observation variances are shared across different state processes; default is all unique
# U.groups is a mx1 vector specifying which u's are shared across different state processes; default is all unique
# Q.groups is a mx1 vector specifying which process variances are shared across different state processes; default is all unique
############# Check that package MASS is installed and load it if it isn't already loaded
###########################################################
if(require(MASS)==FALSE) {
str1 = "The KalmanEM function requires the package MASS and you haven't installed it yet.\n"
str2 = "Install the package before continuing by\n"
str3 = "going to the Packages tab and selecting install (if working from the R GUI).\n"
cat(paste(str1,str2,str3,sep=""))
stop("missing MASS package")
}
############# Set the values for whichPop, Q.groups, U.groups, R.groups and inits if they were not passed in
###########################################################
U.names=U.groups
if(is.numeric(U.groups)==FALSE) U.groups = as.numeric(as.factor(U.groups))
R.names=R.groups
if(is.numeric(R.groups)==FALSE) R.groups = as.numeric(as.factor(R.groups))
Q.names=Q.groups
if(is.numeric(Q.groups)==FALSE) Q.groups = as.numeric(as.factor(Q.groups))
n = dim(y)[1]
default = length(whichPop)==1 & is.na(whichPop[1]); if(default) whichPop = seq(1,n)
m = max(whichPop)
default = length(Q.groups) == 1 & is.na(Q.groups[1]); if(default) Q.groups = seq(1,m)
default = length(U.groups) == 1 & is.na(U.groups[1]); if(default) U.groups = seq(1,m)
default = length(R.groups) == 1 & is.na(R.groups[1]); if(default) R.groups = seq(1,n)
default = length(Uinit)==1 & is.na(Uinit[1]); if(default) Uinit = rep(0,m)
default = length(Qinit)==1 & is.na(Qinit[1]); if(default) Qinit = rep(0.05,m)
default = length(Ainit)==1 & is.na(Ainit[1]); if(default) Ainit = rep(0,n)
default = length(Rinit)==1 & is.na(Rinit[1]); if(default) Rinit = rep(0.05,n)
default = length(x00)==1 & is.na(x00[1]);
tmp = y; tmp[tmp==-99] = NA;
if(default) x00 = rep(mean(tmp,na.rm=T), m) #this use the mean over all the observations for all observation ts
default = length(V00init)==1 & is.na(V00init[1]); if(default) V00init = rep(1,m);
if(length(which(is.na(U.names)==TRUE))) U.names=U.groups
if(length(which(is.na(Q.names)==TRUE))) Q.names=Q.groups
if(length(which(is.na(R.names)==TRUE))) R.names=R.groups
############ Error check to make sure user did not specify an illegal or inconsistent model
###########################################################
if(FailedErrorCheck(y, whichPop, U.groups, R.groups, Q.groups, varcov.Q, varcov.R,Uinit,Qinit,Ainit,Rinit,x00))
stop("Model is illegal or not interally consistent")
############# Set up the model structure: how many state processes and which observation time series go to which state process
###########################################################
n <- dim(y)[1] #the user passes in y with time down the rows, but at top of function this is transposed because this func wants time in the cols
m = max(whichPop) # this is the number of state processes
numYrs <- dim(y)[2]
# Construct the Z matrix based on n, m, and whichPop; this 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
# Set up the grouping structure among different state processes and the different observation time se
U.numGroups <- max(U.groups)
Q.numGroups <- max(Q.groups)
R.numGroups <- max(R.groups)
# Create matrices of 0s and 1s for each, to later be used for multiplication
ZU = matrix(0,m,U.numGroups) # matrix to allow shared growth rates
for(i in 1:U.numGroups) ZU[which(U.groups == i),i] <- 1
ZQ <- matrix(0,m,Q.numGroups) # matrix to allow shared process variances
for(i in 1:Q.numGroups) ZQ[which(Q.groups==i),i] <- 1
ZR <- matrix(0,n,R.numGroups) # matrix to allow shared measurement errs
for(i in 1:R.numGroups) ZR[which(R.groups==i),i] <- 1
# Print out the model structure as a check for the user
if(!silent) {
modelStruc = paste("Model Structure is\n","m: ",m," state process(es)\n","n: ",n," observation time series\n",sep="")
modelStruc = paste(modelStruc,"whichPop: Observation time series assigned to state processes as ",paste(whichPop,collapse=" "),"\n",sep="")
if(m > 1) modelStruc = paste(modelStruc,"U.groups: State process growth rates assigned to groups as ",paste(U.groups,collapse=" "),"\n",sep="")
if(m > 1 & varcov.Q == "unconstrained")
modelStruc = paste(modelStruc,"Q: Process errors have an unconstrained m x m variance-covariance matrix\n",sep="")
if(m > 1 & varcov.Q == "equalvarcov")
modelStruc = paste(modelStruc,"Q: Process errors have an m x m variance-covariance matrix with equal variances and equal covariances\n",sep="")
if(m > 1 & varcov.Q == "diagonal"){
modelStruc = paste(modelStruc,"Q: Process errors are uncorrelated and have a diagonal var-cov matrix.\n",sep="")
modelStruc = paste(modelStruc,"Q.groups: State process variances assigned to groups as ",paste(Q.groups,collapse=" "),"\n",sep="")
}
if(varcov.R == "unconstrained")
modelStruc = paste(modelStruc,"R: Observation errors have an unconstrained n x n variance-covariance matrix\n",sep="")
if(varcov.R == "equalvarcov")
modelStruc = paste(modelStruc,"R: Observation errors have an n x n variance-covariance matrix with equal variances and equal covariances\n",sep="")
if(varcov.R == "diagonal"){
modelStruc = paste(modelStruc,"R: Observation errors are uncorrelated and have a diagonal var-cov matrix.\n",sep="")
modelStruc = paste(modelStruc,"R.groups: Observation variances assigned to groups as ",paste(R.groups,collapse=" "),"\n",sep="")
}
cat(modelStruc)
}
#############Set up some matrices used 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
#############Create record of each variable over all iterations
###########################################################
iter.record <- list(U = array(0,dim=c(m,max.iter)), Q = array(0,dim=c(m,m,max.iter)),
A = array(0,dim=c(n,max.iter)), R = array(0,dim=c(n,n,max.iter)), loglike = array(0, dim=c(1, max.iter)) )
#############Start the EM algorithm which is going to keep updating until
# it converges
###########################################################
#Set the first (or initial) values for the parameters
U= array(Uinit, dim=c(m,1))
Q=makediag(Qinit,nrow=m)
A= array(Ainit,dim=c(n,1))
R=makediag(Rinit,nrow=n) #again we should be estimating this
x00=array(x00,dim=c(m,1))
init.x00 = x00
V00 <- makediag(V00init,nrow=m); #initial variance of initx;
B = makediag(1,m) # Interaction matrix: If estimated it will change from identity, otherwise remain constant
initY1 = 0
for(i in 1:n) {
# do a simple linear regression, estimating the predicted value
y.obs = y[i,]
y.obs[which(y.obs==0)] = NA
initY1[i] = lm(y.obs~seq(1:numYrs))$fitted.values[1]
}
bestLL = -1.0e10 #this is only used if MonteCarloInit == TRUE
if(MonteCarloInit == FALSE) {
# use values passed in by user
numInits = 0
numInitSteps = 0
}
# EW: the first numInits iterations of this loop represent the MCinitialization phase (if MCInit=T),
# the last iteration is the final estimation (if MCInit=F, then goes straight to the last)
for(loop in 1:(numInits + 1)) {
#reset for each loop
loglike.old <- -10000
cvg <- 1+tol
# 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))
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)
init.x00 = array(runif(m,0.75*initY1,1.25*initY1),dim=c(m,1))
x00 = init.x00
}
# Some initial start draws of the params are so awful that loglike will go down. We want to discard these cases.
validDraw = TRUE
# If we've exceeded the initialization phase, use the best params for starting final estimation
if(loop > numInits & MonteCarloInit == TRUE) {
Q = best.Q
R = best.R
U = best.U
x00 = best.x00
}
for(iter in 1:numIter) {
kf <- Kfilter(m,n,numYrs,U,Q,A,R,B,x00,V00,include,y) #call kalman filter to get x(t)|y(1:T) etc estimates
loglike.new = kf$loglik
if(iter>1 && is.finite(loglike.old) == T && is.finite(loglike.new) == T) cvg = loglike.new - loglike.old ## cvg is improvement in likelihood
# If loglike goes down, set validDraw to FALSE
if(cvg<0) validDraw = FALSE
else validDraw = TRUE
# if converged, break out
if(cvg >= 0 & cvg < tol) break
# Store loglike for comparison to new one after parameters are updated
loglike.old = loglike.new
################# M STEP Re-estimate U,Q,A,R,B,X00,V00 via ML given x(t) estimate
################
# Get new A subject to its constraints (update of R will use this)
# See notes by EH for derivation
##############################################################
A.last.iter <- A
A=array(0,dim=c(n,1)); #A will be zero except when there are multiple observation time series for a single state process
# any state process with multiple observation time series will have 1st A set to 0 and rest estimated
tableA <- as.numeric(table(whichPop)) #figure out which state processes have more than 1 observation time series
numSubsWith2Sites = which(tableA > 1)
if(length(numSubsWith2Sites) > 0) { # added by EW 07/19 to catch the case where A not estimated
for(subpops in 1:length(numSubsWith2Sites)) { # loop over all state processes with more than 1 observation time series
this.sub = numSubsWith2Sites[subpops]
this.n = tableA[this.sub]
this.index = which(whichPop==this.sub)
sum1 <- 0
for (i in 1:numYrs) {
include.Y <- include[this.index,this.sub,i] # include is the matrix of 0s/1s for all years, include.Y is the matrix for a single year
dim(include.Y) <- c(this.n,1) #this is just for one state process so nx1
A.if.y.missing = makediag(miss[this.index,i],nrow=this.n)%*%A.last.iter[this.index] #this is going to give me 0 if have val and old A if not
A.if.y.present = array(y[this.index,i],dim=c(this.n,1))-include.Y %*% kf$xtT[this.sub,i]
sum1 <- sum1 + A.if.y.present + A.if.y.missing
} #end for numYrs loop
As = sum1/numYrs
As = As - As[1]
A[this.index] = As
} # end for subpop loop
} # end if
################
# Get new U subject to its constraints (update of Q will use this)
# See notes by EH and Rich Hinrichsen for derivation
################################################################
# if m=1 or m=n this is going to reduce to U <- (kf$xtT[,numYrs]-kf$xtT[,1])/(numYrs-1)
# if some state processes share a u, then we need to take the average across processes sharing a u, taking into account the variance of each process
Qinv = chol2inv(chol(Q)) # this is calculated here because used twice below
Qinv = (Qinv+t(Qinv))/2 #enforce symmetry
numer = t(ZU)%*%Qinv%*%(kf$xtT[,numYrs]-kf$xtT[,1])
denom = 1/takediag((t(ZU)%*%Qinv%*%ZU)*(numYrs-1))
U = c(ZU%*%(numer*denom))
################
# Get new R subject to its constraints (updated 07.09.08 by EH to fix R est bug and by EW to add grouping)
# S&S 4.77 with addition of grouping and diagonal constraint variants
################################################################
sum1 <- 0
for (i in 1:numYrs) {
include.Y <- include[,,i] # include is the matrix of 0s/1s for all years, include.Y is the matrix for a single year
dim(include.Y) <- c(n,m) #stops R from changes the dimensions when m=1
include.Y2 <- array(notmiss[,i],dim=c(n,1))
# this is the updating equation for R
# EH 7.10.08 THIS IS THE PART THAT IS PREVENTING ESTIMATION OF R COVARIANCES WHEN THERE ARE MISSING VALUES
R.if.y.missing <- makediag(miss[,i],nrow=n)%*%R #this is going to give me 0 if have val and R from last iteration if not
err <- y[,i]-(include.Y %*%kf$xtT[,i]+include.Y2*A) #residuals: this will be zero when y[j] is missing]
R.if.y.present <- err%*%t(err) + include.Y %*% kf$VtT[,,i] %*% t(include.Y) #updated R estimate if have y data
sum1 <- sum1 + R.if.y.present + R.if.y.missing
} #end for loop
sum1=(sum1+t(sum1))/2 #enforce symmetry;
R_unconstrained = sum1/numYrs #this provides the estimate of the R matrix with diagonal and non-diagonal elements
#The inv functions don't work well when the diagonal elements get lower than the machine tol, so put a min on those
#this could break if R_unconstrained happens to equal an integer
if(length(R_unconstrained)==1)
R_unconstrained <- max(sqrt(.Machine$double.eps),R_unconstrained)
else
diag(R_unconstrained)[which(takediag(R_unconstrained) 1) #this will get the average to use for the shared covariance value
Roffdiag <- (1/(n*(n-1)))*(sum(R_unconstrained)-sum(takediag(R_unconstrained)))
Roffdiag <- Roffdiag*(matrix(1,n,n)-makediag(1,n))
R <- Rdiag + Roffdiag
} #if R is equal var and equal cov
################
# Get new Q subject to its constraints (updated 07.09.08 by by EW to add grouping)
# This is a little different than S&S Eqn 4.76; theirs is based on the likelihood
# as written in Eqn 4.69 & and Harvey 4.2.19
# We don't have a prior on V00 and what you set that at affects the Q M-step calculation
# to its detriment.
# Instead I'm using the likelihood calculation from Ghahramani and Hinton which is does
# the sum for the Q bit from 2 to T not 1 to T; see notes by EH
# S&S and Harvey would start here
# S00 <- kf$V0T + kf$x0T%*%t(kf$x0T)
# S11 <- kf$VtT[,,1] + (kf$xtT[,1]-U)%*%t(kf$xtT[,1]-U)
# 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
################################################################
S00 = 0
S11 = 0
S10 = 0
for (i in 2:numYrs) {
S00 <- S00 + (kf$VtT[,,i-1] + kf$xtT[,i-1]%*%t(kf$xtT[,i-1]));
S10 <- S10 + (kf$Vtt1T[,,i]+(kf$xtT[,i]-U)%*%t(kf$xtT[,i-1]));
S11 <- S11 + (kf$VtT[,,i] + (kf$xtT[,i]-U)%*%t(kf$xtT[,i]-U));
}
Q_unconstrained = (S11 - (S10 + t(S10)) + S00)/(numYrs-1); #numYrs - 1 since summing 2 to T
# Calculate the interaction matrix, if requested
if(estInteractions==TRUE) B = S10%*%chol2inv(chol(S00))
#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) 1)
Qoffdiag <- (1/(m*(m-1)))*(sum(Q_unconstrained)-sum(takediag(Q_unconstrained)))
Qoffdiag <- Qoffdiag*(matrix(1,m,m)-makediag(1,m))
Q <- Qdiag + Qoffdiag;
} #if equal variances and covariances
################
# EH: Get new x00, V00 cannot be estimated so it is fixed
# S&S Eqn 4.78
################################################################
x00 <- c(kf$x0T)
################
# Keep a record of the iterations; this is to allow user to debug if they are having problems with convergence
################################################################
iter.record$R[,,iter] <- R
iter.record$U[,iter] <- U
iter.record$Q[,,iter] <- Q
iter.record$A[,iter] <- A
iter.record$loglike[iter] <- loglike.new
} # end inner iter loop
# If doing a MCInit, check whether the likelihood is the best observed
# Only use bootstrap param draws where loglike did not go down during numInitSteps
if(validDraw == TRUE & MonteCarloInit == TRUE & loglike.new > bestLL) {
# update the best initial parameter estimates
best.Q = iter.record$Q[,,1]
best.R = iter.record$R[,,1]
best.U = iter.record$U[,1]
best.x00 = init.x00
bestLL = loglike.new
}
} # end iter loop that is running until the EM algorithm converges
# Consolidate arrays for output
iter.record$R <- iter.record$R[,,1:(iter-1)]
iter.record$U <- iter.record$U[,1:(iter-1)]
iter.record$Q <- iter.record$Q[,,1:(iter-1)]
iter.record$A <- iter.record$A[,1:(iter-1)]
iter.record$loglike <- iter.record$loglike[1:(iter-1)]
if(!silent){
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])))
}
names(U) = U.names
if(is.matrix(Q)==F) Q = as.matrix(Q) # this line added to make sure that a scalar isn't passed to rmvnorm - 01/20/09
rownames(Q) = Q.names
colnames(Q) = Q.names
rownames(R) = R.names
colnames(R) = R.names
states = kf$xtT
return(list(states=states,CI_states.lower = (kf$xtT - 1.96*se), CI_states.upper = (kf$xtT + 1.96*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.new,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, data=t(y), include=include))
} #End of the KalmanEM function
Kfilter = function(m,n,numYrs,U,Q,A,R,Phi,x00,V00,include,y) {
# EH: This is running a Kalman-Raush filter (forward + backward) to get (written univariate to make it look cleaner):
# xtT: E(x(t) | y(1:T)) , VtT: Var(x(t)*x(t) | y(1:T)), Vtt1T: Var(x(t)*x(t-1) | y(1:T))
# EJW: This code is a hybrid from Shumway and Stoffer (2006) and Matlab
# code originally written by Eli Holmes for multi-state state-space
# estimation. The code originally used solve() for matrix inverses -
# EW changed this to cholesky decomp though for speed.
# ** All eqn refs are to 2nd ed of Shumway & Stoffer (2006): Time Series Analysis and Its Applications
#initialize - these are for the forward, Kalman, filter
# for notation purposes, 't' represents current point in time, 'T' represents the length of the series
Vtt <- array(0,dim=c(m,m,numYrs)) # Analagous to S&S Ptt, E[Vtt|y(1:t)]
Vtt1 <- array(0,dim=c(m,m,numYrs)) # Analagous to S&S Ptt1, E[Vtt1|y(1:t)]
xtt <- array(0,dim=c(m,numYrs)) # E[x(t) | Y(t)]
xtt1 <- array(0,dim=c(m,numYrs)) # E[x(t) | Y(t-1)]
innov <- array(0,dim=c(n,numYrs)) # these are innovations, parentheses in 6.21
vt <- array(0,dim=c(n,numYrs)) # used for likelihood, vt equivalent to epsilon, eqn 6.62
Ft <- array(0,dim=c(n,n,numYrs)) # used for likelihood, Ft equivalent to sigma matrix eqn 6.62
# these are for the backwards, Rausch, filter
VtT <- array(0,dim=c(m,m,numYrs)) # E[Vtt|y(1:T)]
J <- array(0,dim=c(m,m,numYrs)) # see eqn 6.49
Vtt1T <- array(0,dim=c(m,m,numYrs)) # E[Vtt1|y(1:T)]
xtT <- array(0,dim=c(m,numYrs)) # E[x | y(1:T)]
Kt <- array(0, dim=c(m,n,numYrs)) # 3D matrix of Kalman gain, EW added 11/14/08
#forward pass gets you E[x(t) given y(1:t)]
xtt1[,1] <- Phi%*%c(x00) + U # eqn 6.19
Vtt1[,,1] <- Phi%*%V00%*%t(Phi) + Q # eqn 6.20
include.Y <- include[,,1] # used for missing data
dim(include.Y) <- c(n,m) #stops R from changes the dimensions when m=1
include.Y2 <- (include.Y%*%array(1,dim=c(m,1))) # EH: create include.Y for A that is nx1
siginv = include.Y%*%Vtt1[,,1]%*%t(include.Y)+R # bracketed piece of eqn 6.23
siginv=chol2inv(chol(siginv)) # now siginv is sig[[i]]^{-1}
siginv = (t(siginv)+siginv)/2 #Vtt1 happens to be symmetric since it is V00+Q; although in general E(xt t(xt1)) is not symmetric
Kt[,,1] <- Vtt1[,,1]%*%t(include.Y) %*% siginv; #broke siginv to impose symmetry, eqn 6.23
innov[,1] <- y[,1] - (include.Y%*%xtt1[,1] + include.Y2*A)
xtt[,1] <- xtt1[,1] + Kt[,,1]%*%innov[,1]; # eqn 6.21
Vtt[,,1] <- Vtt1[,,1]-Kt[,,1]%*%include.Y%*%Vtt1[,,1]; # eqn 6.22, detail after 6.28
Vtt[,,1] <- (Vtt[,,1]+t(Vtt[,,1]))/2; #Vtt must be symmetric
# the missing values will contribute 0.0 for this calc
vt[,1] <- y[,1]-(include.Y%*%xtt1[,1]+include.Y2*A) #3.3.4; need to hold on to this for loglik calc
Ft[,,1] <- include.Y%*%Vtt1[,,1]%*%t(include.Y)+R #need to hold on to this for loglik calc
Ft[,,1] <- (Ft[,,1]+t(Ft[,,1]))/2; #to ensure its symetric
for (t in 2:numYrs) {
xtt1[,t] <- Phi%*%xtt[,t-1] + U #xtt1 denotes x_t^(t-1), eqn 6.19
Vtt1[,,t] <- Phi%*%Vtt[,,t-1]%*%t(Phi) + Q # eqn 6.20
Vtt1[,,t] <- (Vtt1[,,t]+t(Vtt1[,,t]))/2 #in general Vtt1 is not symmetric but here it is since Vtt and Q are
include.Y <- include[,,t]
dim(include.Y) <- c(n,m) #stops R from changes the dimensions when m=1
include.Y2 <- (include.Y%*%array(1,dim=c(m,1))) # I'm futzing with the include.Y for A to make it nx1
Kt[,,t] <- Vtt1[,,t]%*%t(include.Y) %*% chol2inv(chol(include.Y%*%Vtt1[,,t]%*%t(include.Y)+R)) # eqn 6.23
innov[,t] <- y[,t] - (include.Y%*%xtt1[,t] + include.Y2*A)
xtt[,t] <- xtt1[,t] + Kt[,,t]%*%innov[,t] # eqn 6.21
Vtt[,,t] <- Vtt1[,,t]-Kt[,,t]%*%include.Y%*%Vtt1[,,t] # eqn 6.22, detail after 6.28
Vtt[,,t] <- (Vtt[,,t]+t(Vtt[,,t]))/2 #to ensure its symetric
# missing values will contribute nothing
vt[,t] <- y[,t]-(include.Y%*%xtt1[,t]+include.Y2*A) #need to hold on to this for loglik calc
Ft[,,t] <- include.Y%*%Vtt1[,,t]%*%t(include.Y)+R #need to hold on to this for loglik calc
Ft[,,t] <- (Ft[,,t]+t(Ft[,,t]))/2 #to ensure its symetric
}
KT <- Kt[,,t];
#indexing is 0 to T for the backwards (Rausch) recursions
#backward pass gets you E[x(t)|y(1:T)] from E[x(t)|y(1:t)]
xtT[,numYrs] <- xtt[,numYrs]
VtT[,,numYrs] <- Vtt[,,numYrs]
s <- seq(numYrs,2)
for(i in 1:(numYrs-1)) {
yr <- s[i] #equivalent to T:-1:0
Vinv <- chol2inv(chol(Vtt1[,,yr]))
Vinv <- (Vinv + t(Vinv))/2 #to enforce symmetry after chol2inv call
J[,,yr-1] <- Vtt[,,yr-1]%*%t(Phi)%*%Vinv # eqn 6.49
xtT[,yr-1] <- xtt[,yr-1] + J[,,yr-1]%*%(xtT[,yr]-xtt1[,yr]) # eqn 6.47
VtT[,,yr-1] <- Vtt[,,yr-1] + J[,,yr-1]%*%(VtT[,,yr]-Vtt1[,,yr])%*%t(J[,,yr-1]) # eqn 6.48
VtT[,,yr-1] <- (VtT[,,yr-1]+t(VtT[,,yr-1]))/2 #VtT is symmetric
}
Vinv <- chol2inv(chol(Vtt1[,,1]))
Vinv <- (Vinv + t(Vinv))/2 #to enforce symmetry after chol2inv call
J0 <- V00%*%t(Phi)%*%Vinv # eqn 6.49
x0T <- x00 + J0%*%(xtT[,1]-xtt1[,1]); # eqn 6.47
V0T <- V00 + J0%*%(VtT[,,1]-Vtt1[,,1])*t(J0) # eqn 6.48
V0T <- (V0T+t(V0T))/2;
#run another backward recursion to get E[x(t)x(t-1)|y(T)]
include.Y <- include[,,numYrs]
dim(include.Y) <- c(n,m) #stops R from changes the dimensions when m=1
Vtt1T[,,numYrs] <- (makediag(1,m) - KT%*%include.Y)%*%Vtt[,,numYrs-1] #this is Var(x(T)x(T-1)|y(T))
s <- seq(numYrs,3)
for (i in 1:(numYrs-2)) {
yr <- s[i]
Vtt1T[,,yr-1] <- Vtt[,,yr-1]%*%t(J[,,yr-2]) + J[,,yr-1]%*%(Vtt1T[,,yr]-Vtt[,,yr-1])%*%t(J[,,yr-2])
}
Vtt1T[,,1] <- Vtt[,,1]%*%t(J0) + J[,,1]%*%(Vtt1T[,,2]-Vtt[,,1])%*%t(J0)
#Calculate log likelihood, see eqn 6.62
loglik <- -n*(numYrs/2)*log(2*pi)
for (i in 1:numYrs) {
if(length(Ft[,,i])==1) detFt <- Ft[,,i] else detFt <- det(Ft[,,i])
Ftinv <- chol2inv(chol(Ft[,,i]))
Ftinv <- (Ftinv +t(Ftinv))/2 #enforce symmetry; Ft is symmetric
loglik <- loglik - (1/2)%*%t(vt[,i]) %*% Ftinv %*% vt[,i] - (1/2)*log(detFt);
}
return(list(xtT = xtT, VtT = VtT, Vtt1T = Vtt1T, x0T = x0T, V0T = V0T, loglik = loglik, Vtt = Vtt, Vtt1 = Vtt1, J=J, Kt=Kt, xtt1 = xtt1, Innov=innov, Sigma=Ft))
}
FailedErrorCheck = function(y, whichPop, U.groups, R.groups, Q.groups, varcov.Q, varcov.R,Uinit,Qinit,Ainit,Rinit,x00) {
############# Error check to make sure user isn't trying to set up an illegal or inconsistent model structure
#################### Check whichPop
n = dim(y)[1] #each row is obs ts
if(length(whichPop) != n) {
cat("Error: whichPop is not in the correct form; it should be 1xn or nx1; \n"); return(TRUE) }
m = max(whichPop) # this is the number of state processes
if( FALSE %in% ( (1:m) %in% whichPop ) ) {
cat("Error: Something is wrong with whichPop. You need an observation time series for each state process. \n"); return(TRUE) }
#################### Check varcov.Q
if( (varcov.Q %in% list("unconstrained", "diagonal", "equalvarcov"))==FALSE ) {
cat("Error: options for varcov.Q are scalar, unconstrained, diagonal or equalvarcov (passed as text in quotes) \n" ); return(TRUE) }
#################### Check varcov.R
if( (varcov.R %in% list("unconstrained", "diagonal", "equalvarcov"))==FALSE ) {
cat("Error: options for varcov.R are unconstrained, diagonal or equalvarcov (passed as text in quotes) \n" ); return(TRUE) }
if( (-99 %in% y) & varcov.R != "diagonal" ) {
cat("Sorry! The current code cannot estimate observation error covariances when there are missing data; set varcov.R to diagonal \n"); return(TRUE) }
#################### Check Q.groups
if(length(Q.groups) != m) {
cat("Error: There is a mismatch between the number of state processes and length of Q.groups \n"); return(TRUE) }
not.all.m.in.groups = FALSE %in% ( (1:m) %in% Q.groups )
if( not.all.m.in.groups & varcov.Q != "diagonal" ) {
cat("Sorry! The code cannot apply groupings of Qs except to diagonal Q matrices \n"); return(TRUE) }
# Check that each state process is assigned to a Q group
if(FALSE %in% ((1:max(Q.groups)) %in% Q.groups) ) {
cat(paste("Error: something is wrong with Q.groups You have ",max(Q.groups)," Q groups but not all appear in Q.groups \n",sep="")); return(TRUE) }
#################### Check U.groups
if(length(U.groups) != m) {
cat(paste("Error: U.groups should be m x 1 or 1 x m","Model is m = ",m," as currently set up.\n",sep="")); return(TRUE) }
if(FALSE %in% ((1:max(U.groups)) %in% U.groups) ) {
cat(paste("Error: something is wrong with U.groups. You have ",max(U.groups)," U groups but not all appear in U.groups. \n",sep="")); return(TRUE) }
#################### Check R.groups
if(length(R.groups) != n) {
cat("Error: There is a mismatch between the number of measurement time series and length of R.groups. \n"); return(TRUE) }
not.all.n.in.groups = FALSE %in% ( (1:n) %in% R.groups )
if( not.all.n.in.groups & varcov.R != "diagonal" ) {
cat("Sorry! The code cannot apply R.groups except to diagonal R matrices; set varcov.R to diagonal. \n"); return(TRUE) }
if(FALSE %in% ((1:max(R.groups) ) %in% R.groups) ) {
cat(paste("Error: something is wrong with R.groups. You have ",max(R.groups)," R groups but not all appear in R.groups.\n",sep="")); return(TRUE) }
############# Check that the initial values are correctly sized
###########################################################
if(length(Uinit) != m) {
cat("Error: Uinit should be m x 1 or 1 x m; Either your Uinit is wrong or your flags are setting m to something you are not expecting.\n"); return(TRUE) }
if(length(Qinit) != m) {
cat("Error: Qinit should be m x 1 or 1 x m; Either your Qinit is wrong or your flags are setting m to something you are not expecting.\n"); return(TRUE) }
if(length(Ainit) != n) {
cat("Error:L Ainit should be n x 1 or 1 x n; Your Ainit is the wrong size.\n"); return(TRUE) }
if(length(Rinit) != n) {
cat("Error: Rinit should be n x 1 or 1 x n; Your Rinit is the wrong size.\n"); return(TRUE) }
if(length(x00) != m) {
cat("Error: x00 should be m x 1 or 1 x m; Either your x00 is wrong or your flags are setting m to something you are not expecting.\n"); return(TRUE) }
return(FALSE)
}
takediag = function(x)
{
if(length(x)==1) return(x)
else return(diag(x))
}
makediag = function(x,nrow=NA)
{
if(length(x)==1)
{
if(is.na(nrow)) nrow=1
return(diag(c(x),nrow))
}
if((is.matrix(x) | is.array(x)))
if(!(dim(x)[1]==1 | dim(x)[2]==1)) stop("Error in call to makediag; x is not vector")
if(is.na(nrow)) nrow=length(x)
return(diag(c(x),nrow))
}
#######################################################################################################
# These functions are all for bootstrapping
#######################################################################################################
bootstrapMSSM = function(nBoot, model, maxIter = 5000, tol=0.01, onlyAIC = TRUE, silent=FALSE) {
# This function updated by EW and EH
# all equation numbers refer to Chapter 6, Shumway & Stoffer
# onlyAIC is a debug feature to allows all the boot params and boot data iterates to be returned
# silent is a flag to indicate whether the progress bar should be printed
###### Check that package mvtnorm is installed and load it if it isn't already loaded
if(require(mvtnorm)==FALSE) {
str1 = "The bootstrapping functions require the package mvtnorm and you haven't installed it yet.\n"
str2 = "Install the package before continuing by\n"
str3 = "going to the Packages tab and selecting install (if working from the R GUI).\n"
cat(paste(str1,str2,str3,sep=""))
stop("the mvtnorm package is not installed")
}
errors = t(stdInnov(model$Sigma, model$Innov)) # standardize innovations
numStates = model$m # number of subpopulations
numYears = model$numYrs # length of time series
numSites = model$n # number of time series
# Calculate the sqrt of sigma matrices, so they don't have to be computed 5000+ times
sigma.Sqrt = array(0,dim=c(numSites, numSites, numYears))
BKS = array(0,dim=c(numStates, numSites, numYears))
for(i in 1:numYears) {
eig = eigen(model$Sigma[,,i]) # calculate inv sqrt(sigma[1])
# EW added as.matrix 01/26/09
sigma.Sqrt[,,i] = eig$vectors %*% makediag((sqrt(eig$values))) %*% t(eig$vectors)
BKS[,,i] = model$B %*% model$Kt[,,i] %*% sigma.Sqrt[,,i] # this is m x n
}
U=NA; Q=NA; R=NA; A=NA; B=NA; states=NA; boot.data=NA;
# these are the arrays for output
if(onlyAIC == FALSE) {
boot.data = array(0,dim=c(dim(as.matrix(model$data)),nBoot))
U = array(0,dim=c(dim(as.matrix(model$U)),nBoot))
Q = array(0,dim=c(dim(as.matrix(model$Q)),nBoot))
R = array(0,dim=c(dim(as.matrix(model$R)),nBoot))
A = array(0,dim=c(dim(as.matrix(model$A)),nBoot))
B = array(0,dim=c(dim(as.matrix(model$B)),nBoot))
states = array(0,dim=c(dim(as.matrix(model$states)),nBoot))
}
logL.thetaStar = 0
logL = 0;
drawProgressBar = FALSE #If not time library, no prog bar
if(require(time)==TRUE) { #then we can draw a progress bar
prev = progressBar()
drawProgressBar = TRUE
}
for(i in 1:nBoot) {
newData = makeNewData(errors, model, model$Sigma, sigma.Sqrt, BKS) # make new data
# Take the ML output and make sure it's in the right form to be passed back in
Rstart = diag(model$R)
Qstart = diag(model$Q)
# Step 2: estimate the MLEs, theta*
boot.model = KalmanEM(newData, varcov.Q=model$varcov.Q, varcov.R=model$varcov.R, whichPop=model$whichPop,U.groups=model$U.groups, Q.groups=model$Q.groups,R.groups=model$R.groups,Uinit=model$U,Ainit=model$A,Qinit=Qstart,Rinit=Rstart,tol=tol,silent=T,max.iter=maxIter)
if(onlyAIC == FALSE) {
# Step 3: store the parameter estimates from the bootstraped data fitting
boot.data[,,i] = as.matrix(newData)
U[,,i] = as.matrix(boot.model$U)
Q[,,i] = as.matrix(boot.model$Q)
R[,,i] = as.matrix(boot.model$R)
A[,,i] = as.matrix(boot.model$A)
B[,,i] = as.matrix(boot.model$B)
states[,,i] = as.matrix(boot.model$states)
}
logL[i] = Kfilter(numStates, numSites, numYears, boot.model$U, boot.model$Q, boot.model$A, boot.model$R, boot.model$B, boot.model$x00, boot.model$V00, boot.model$include, t(model$data))$loglik
if(drawProgressBar) prev <- progressBar(i/nBoot,prev)
}
AICb = -4*(sum(logL))/nBoot + 2*model$loglike
return(list(AICb = AICb, logL.star=logL, U = U, Q = Q, R = R, A = A, B = B, states = states, boot.data=boot.data))
}
makeNewData = function(Errors, model, Sigma, sigmaSqrt, B.Kt.sqrtSigma) {
# This function updated by EW and EH Dec 02, 2008
# Errors = std innovations
origData = model$data
newData = matrix(-99, model$numYrs, model$n)
newStates = matrix(-99, model$numYrs+1, model$m) # States = years x subpops
numMissingVals = length(which(Errors==0))
if(numMissingVals == 0) {
# Then the bootstrap algorithm proceeds
# Stoffer & Wall suggest not sampling from innovations 1-3 (not much data)
minIndx = ifelse(model$numYrs > 5, 4, 1)
samp = sample(seq(minIndx, model$numYrs), size = model$numYrs, replace = T)
# EW added as.matrix() on 01/26/09
e = as.matrix(Errors[samp,])
# newStates is a[] in the writeup by EH
newStates[1,] = model$x00
#newStates[2,] = model$B %*% model$x00 + model$U + model$B %*% model$Kt[,,1] %*% sigmaSqrt[,,1] %*% e[1,]
newStates[2,] = model$B %*% model$x00 + model$U + B.Kt.sqrtSigma[,,1] %*% e[1,]
newData[1,] = model$x00[model$whichPop] + model$A + sigmaSqrt[,,1] %*% e[1,] # eqn 1.7 in writeup by EH
for(i in 2:model$numYrs) { # Deal with years 2:T, use sigma, Kt, B, U, A
newStates[i+1,] = model$B %*% newStates[i,] + model$U + B.Kt.sqrtSigma[,,i] %*% e[i,]
newData[i,] = newStates[i,][model$whichPop] + model$A + sigmaSqrt[,,i] %*% e[i,]
}
newData[which(e == 0)]=-99 #newData is going to be used in the KalmanEM function which expects -99 for missing values
# return the new data
}
if(numMissingVals != 0) {
newStates[1,] = model$x00 # t = 0
# create a matrix of RVs for observation error
obs.error = rmvnorm(model$numYrs, mean = rep(0, nrow(as.matrix(model$R))), sigma = model$R, method="chol")
pro.error = rmvnorm(model$numYrs, mean = rep(0, nrow(model$Q)), sigma = model$Q, method="chol")
for(i in 2:(model$numYrs+1)) {
newStates[i,] = model$B %*% newStates[i-1,] + model$U + pro.error[i-1,]
newData[i-1,] = model$Z %*% newStates[i,] + model$A + obs.error[i-1,]
}
newData[which(Errors==0)]=-99
}
return(newData)
}
stdInnov = function(SIGMA, INNOV) {
# This function added by EW Nov 3, 2008
# This function is only called once by bootstrapMSSM()
# SIGMA is covariance matrix, E are original innovations
numYrs = dim(INNOV)[2]
numSites = dim(INNOV)[1]
SI = matrix(0, numSites, numYrs)
for(i in 1:numYrs) {
a = SIGMA[,,i]
a.eig = eigen(a)
# EW added as.matrix() on 01/26/09 for single trajectory model
a.sqrt = a.eig$vectors %*% makediag((sqrt(a.eig$values))) %*% solve(a.eig$vectors)
SI[,i] = chol2inv(chol(a.sqrt)) %*% INNOV[,i] # eqn 1, p359 S&S
}
SI[which(INNOV == 0)]=0 # make sure things that should be 0 are 0
return(SI)
}