KalmanEM = function(y, whichPop=NA, Uinit=NA, Qinit=NA, Ainit=NA, Rinit=NA, x00=NA, V00init=NA, max.iter=100, varcov.Q="unconstrained", varcov.R="diagonal", U.groups=NA, Q.groups=NA, R.groups=NA, U.bounds=c(-1,1), logQ.bounds = c(log(0.0001),log(1.0)), logR.bounds = c(log(0.0001),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
############## 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
###########################################################
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);
############ 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
y[which(y == -99)] <- 0 # Clear out -99s in data; include has been set up so that these values will be ignored
miss <- ifelse(abs(y)>0,0,1) # 0=observed, 1=missing
notmiss <- ifelse(abs(y)>0,1,0) # 1=observed, 0=missing
#############Set up some counters and containers that are used during each iteration
###########################################################
LL <- NA;
converged <- 0
previous_loglik <- -Inf
num_iter <- 0
iter <- 0
loglike.old <- -10000
cvg <- 1+tol
# 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
# EW: the first numInits iterations of this loop represent the initialization phase,
# the last iteration is the final estimation
if(MonteCarloInit == FALSE) {
# use values passed in by user
numInits = 0
numInitSteps = 0
}
for(loop in 1:(numInits + 1)) {
numIter = ifelse(loop <= numInits,numInitSteps,max.iter) # this determines # of iterations, depends on phase
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
}
validDraw = TRUE
# If we've exceeded the initialization phase, do the 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(cvg<0) {
#cat(paste("Warning: at iter=",iter," likelihood decreased; loglike.new - loglike.old = ", cvg,"\n", sep=""))
validDraw = FALSE
}
#if(cvg >= 0 & cvg < tol) break
if(cvg >= 0 & cvg < tol & loop > numInits) break # only break if on the last run
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
# After this iteration, check whether the likelihood is the best observed
if(validDraw == TRUE & loglike.new > bestLL) {
# update the best 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])))
}
return(list(states=kf$xtT,CI_states.lower = (kf$xtT - 1.96*se), CI_states.upper = (kf$xtT + 1.96*se), U=U,Q=Q,A=A,R=R,B=B,x00=x00,V00=V00,V0T=kf$V0T,iter=iter,loglike=loglike.new,iter.record=iter.record))
} #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)]
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)]
#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 <- Vtt1[,,1]%*%t(include.Y) %*% siginv; #broke siginv to impose symmetry, eqn 6.23
xtt[,1] <- xtt1[,1] + Kt%*%(y[,1] - (include.Y%*%xtt1[,1] + include.Y2*A)); # eqn 6.21
Vtt[,,1] <- Vtt1[,,1]-Kt%*%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 <- Vtt1[,,t]%*%t(include.Y) %*% chol2inv(chol(include.Y%*%Vtt1[,,t]%*%t(include.Y)+R)) # eqn 6.23
xtt[,t] <- xtt1[,t] + Kt%*%(y[,t] - (include.Y%*%xtt1[,t] + include.Y2*A)) # eqn 6.21
Vtt[,,t] <- Vtt1[,,t]-Kt%*%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;
#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))
}
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))
}