#Utility function to generate data GenerateData = function(nTime,npad,np,noStates,noObs,r,psi,lambda,theta, Phi_SS,eta_intercept,G_SS,a0,P0,zt, fileObs=NULL, fileState=NULL,shocks=NULL){ yall <- matrix(0,nrow = nTime*np, ncol = noObs) etaall <- matrix(0,nrow = nTime*np, ncol = noStates) if (is.null(shocks)){ shocks = matrix(0,np*nTime,noStates) } for (p in 1:np){ # Set up matrix for state variables. etaC <- matrix(0, nrow = noStates, ncol = nTime + npad) etaC[,1] <- a0 + chol(P0)%*%rnorm(noStates) #Process noises zeta <- mvtnorm::rmvnorm(nTime+npad, mean = rep(0,noStates), sigma = psi) # Measurement errors. epsilon <- rmvnorm(nTime+npad, mean = rep(0,noObs), sigma = theta) # Generate latent state values for (i in 2:(nTime+npad)){ etaC[ ,i] <- eta_intercept + Phi_SS %*% etaC[ ,i-1] + G_SS%*%zt[(i+(p-1)*nTime-1),] + zeta[i, ] + shocks[(i+(p-1)*nTime),] }#End of loop over nTime eta <- t(etaC[,(npad+1):(npad+nTime)]) # generate observed series y <- matrix(0, nrow = nTime, ncol = noObs) for (i in 1:nrow(y)){ y[i, ] <- lambda %*% eta[i, ] + epsilon[i, ] } yall[(1+(p-1)*nTime):(p*nTime),] = y etaall[(1+(p-1)*nTime):(p*nTime),] = eta } #End of p loop yall = cbind(rep(1:np,each=nTime),rep(1:nTime,np),yall[,1:noObs],zt) yall = data.frame(yall) colnames(yall) = c("ID","Time",paste0("V",1:noObs),paste0("z",1:noInput)) etaall = cbind(rep(1:np,each=nTime),rep(1:nTime,np),etaall) colnames(etaall) = c("ID","Time",paste0("LV",1:noStates)) etaall = data.frame(etaall) mytheme <- theme_classic() %+replace% theme(axis.title.x = element_blank(), axis.title.y = element_text(face="bold",angle=90)) p1 <- ggplot(data = etaall, aes(x = Time, y = LV1, group = ID, colour = ID)) + mytheme + labs(title = "Simulated LV", y = "LV 1", x = "nTimeime") + geom_line(size=1) + theme(legend.position="none") p2 <- ggplot(data = etaall, aes(x = Time, y = LV2, group = ID, colour = ID)) + mytheme + labs(title = "Simulated LV", y = "LV 2", x = "nTimeime") + geom_line(size=1) + theme(legend.position="none") p3 <- ggplot(data = yall, aes(x = Time, y = V1, group = ID, colour = ID)) + mytheme + labs(title = "Simulated data", y = "Observed y1", x = "nTimeime") + geom_line(size=1) + theme(legend.position="none") p4 <- ggplot(data = yall, aes(x = Time, y = V2, group = ID, colour = ID)) + mytheme + labs(title = "Simulated data", y = "Observed y2", x = "Time") + geom_line(size=1) + theme(legend.position="none") grid.arrange(p1, p2, p3, p4, ncol = 2) if (!is.null(fileObs)){ write.table(yall,fileObs,row.names=FALSE,col.names=FALSE,sep=",")} if (!is.null(fileState)){ write.table(etaall,fileState,row.names=FALSE,col.names=FALSE,sep=",")} return(list(yall = yall, etaall = etaall)) } #Function to run the controller runController = function(noStates, noObs, noInput, noSubj, noTime, data, idName,timeIndexName,obsNames, inputNames, dynamics, meas, noise, initial, parameterValues, eta_target, Q, Qf, R, h, G_SS2){ I_noStates = diag(rep(1,noStates)) X_NoControl = matrix(NA, ncol = noStates, nrow = noSubj*noTime) #Place holder for latent variables newData = data #Make a copy of the original data passed in for (i in 1:noSubj){ oldInputNames = paste0("old",inputNames) temp.dat = data[(1+(i-1)*noTime):(i*noTime),] temp.dat = cbind(temp.dat,rep(NA,noTime)) colnames(temp.dat)[length(colnames(temp.dat))] = oldInputNames #Set up dynr data and model thedat.original = dynr.data(temp.dat, id=idName, time=timeIndexName, observed=obsNames, covariates=inputNames) model.Nocontrol <- dynr.model(dynamics=dynamics, measurement=meas,noise=noise, initial=initial, data=thedat.original,outfile="temp.c") coef(model.Nocontrol) <- parameterValues res.Nocontrol = dynr.cook(model.Nocontrol, optimization_flag = FALSE,hessian_flag = FALSE, verbose = FALSE, debug_flag=TRUE) X_Originali = res.Nocontrol@eta_filtered X_NoControl[(1+(i-1)*noTime):(i*noTime),] = as.matrix(t(X_Originali),byrow=TRUE,ncol=noStates) temp.dat[,oldInputNames] = temp.dat[,inputNames] L <-list(matrix(0,ncol=noStates,nrow=noStates),h) g <- list(matrix(0,ncol=1,nrow=noStates),h) zstar <- matrix(0,nrow=noInput,ncol=1) L[[h]] <- Q_f # 7a g[[h]] <- -Q_f %*% eta_target[1,] #7b for (t in 1:(noTime-1)){ #If you want the effects of the control inSubjut to affect the dependent #variables in real time, then you need to put the new data #reflecting the effects of the control inSubjut somewhere here, and #re-cook the model to get the state estimates similar to what is #done in the code below that: # Start of "real-time" state estimation # thedat.test = dynr.data(temp.dat[1:(min(noTime,t+h)+(i-1)*noTime),], id=idName, time=timeIndexName, # observed=obsNames,covariates=oldInputNames) # model.test <- dynr.model(dynamics=dynamics, measurement=meas, # noise=noise, initial=initial, data=thedat.test, # outfile="temp0.c") # coef(model.test) <- coef(res) # res.test = dynr.cook(model.test,optimization_flag = FALSE, # hessian_flag = FALSE, # verbose = FALSE, debug_flag=TRUE) # X_OriginalitPlus = res.test@eta_filtered #End of "real-time" state estimation for (tau in (h-1):1) { inSubjutMultiplier <- solve( I_noStates + (L[[tau+1]] %*% G_SS2 %*% solve(R)%*% t(G_SS2))) L[[tau]] <- t(Phi_SS) %*% inSubjutMultiplier%*% L[[tau+1]] %*%Phi_SS + Q g[[tau]] <- t(Phi_SS) %*% inSubjutMultiplier%*% g[[tau+1]] - Q %*% eta_target[min(t+tau,noTime),] if (tau==1){ X_OriginalitPlus = matrix(X_Originali[,t+tau],ncol=1) zstar <- -solve(R) %*% t(G_SS2) %*% inSubjutMultiplier %*% (L[[tau+1]] %*% Phi_SS %*% X_OriginalitPlus + g[[tau+1]]) }} temp.dat[t,inputNames] = matrix(zstar,nrow=1) #Usually this should be [t], set to [t+1] because of #dynr setting for covariates in the dynamic model. temp.dat[t+1,inputNames] = zstar } #noTime-loop newData[(1+(i-1)*noTime):(i*noTime),] = temp.dat[,!colnames(temp.dat) %in% oldInputNames] } #One more pass to get everyone's smoothed state estimates #with zstar thedat.Controlled = dynr.data(newData, id=idName, time=timeIndexName, observed=obsNames,covariates=inputNames) model.Controlled <- dynr.model(dynamics=dynamics, measurement=meas, noise=noise, initial=initial, data=thedat.Controlled, outfile="temp.c") coef(model.Controlled) <- parameterValues res.Controlled <- dynr.cook(model.Controlled,optimization_flag = FALSE, hessian_flag = FALSE, verbose = FALSE, debug_flag=TRUE) X_Controlled <- as.matrix(t(res.Controlled@eta_filtered), ncol=noStates,byrow=TRUE) #Extract filtered estimates #Extract zstar zControlled = as.matrix(newData[,inputNames],col=1) return(list(optimalInput = zControlled, filteredStateswithControl = X_Controlled, filteredStateswithoutControl = X_NoControl, newData = newData)) } #End of controller function # ---- Cost and Benefit Computations ---- # X = Latent variable vector (number of LVs x number of time points) # u = Control inputs (number of inputs x number of time points) # Q = Penalty cost matrix for state deviations # R = Penalty cost matrix for input administration # XLast = Latent Variable vector at terminal time point (number of LVs x 1) # QLast = Penalty cost matrix for state deviations at terminal point # Target = Target state values (number of LVs x number of time points) # Xbase = Baseline latent variable values (usually those from without any control inputs) # ubase = Baseline inputs # Returned output: # A list with elements: # cost = quadratic cost associated with input administration # benefit = proportion of reduction in total state cost associated with (X, XLast) # compared to that of Xbase; # positive value means there is a reduction in this cost of (X,Xlast) relative to Xbase; # negative value means (X,Xlast) is actually doing worse. CostBenefit = function(X, u, Q, R, XLast=NULL, QLast=Q,Target,Xbase=NULL, ubase=NULL){ benefit = 0; cost = 0; nTime = dim(X)[2]; benefitBase=0; costBase=0 for (t in 1:(nTime-1)){ benefit = benefit + as.matrix(t(X[,t] - Target[,t]),nrow=1)%*% Q %*% as.matrix((X[,t]-Target[,t]),ncol=1) if (!is.null(XLast)){ benefit = benefit + as.matrix(t(X[,nTime] - Target[,nTime])%*% QLast %*% as.matrix(X[,nTime]-Target[,nTime]),ncol=1) } if (!is.null(Xbase)){ benefitBase = benefitBase + as.matrix(t(Xbase[,t] - Target[,t]),nrow=1)%*% Q %*% as.matrix((Xbase[,t]-Target[,t]),ncol=1) } cost = cost + as.matrix(t(u[,t]),nrow=1) %*% R %*% as.matrix(u[,t],ncol=1) if (!is.null(ubase)){ costBase = costBase + as.matrix(t(ubase[,t]),nrow=1) %*% R %*% as.matrix(ubase[,t],ncol=1) } } return(list(cost = (cost-costBase)/max(1,costBase), benefit = (benefitBase-benefit)/benefitBase)) }