Lorenz system estimated different ways with ctsem

lorenzfit.R

#get ctsem
install.packages("devtools")
library(devtools)
install_github("cdriveraus/ctsem")

#load lorenz data
ldat=read.csv(file='ldat.csv')
ldat=ldat[seq(1,nrow(ldat),5),] #subsample


#need time and id columns in data
ldat$time <- ldat$X.time
ldat$id <- 1 #subject id

require(ctsem)

#####ctsem specification using time varying par

#dynamic error matrix
DIFFUSION=diag(1e-5,4)
DIFFUSION[4,4] = 'tadynerrorsd'

#model spec (some unnecessary matrices fixed to arbitrary values)
cm <- ctModel(n.latent=4, n.TDpred = 0,n.manifest=1,type='stanct',
  manifestNames = 'X_obs.value',
  DRIFT=diag(-1,4), #not used
  DIFFUSION=DIFFUSION,
  LAMBDA=matrix(c(1,0,0,0),1), #effect of latent states on measurements
  CINT=matrix(0,4,1), #not used
  PARS = matrix(c('tadyn','p','B', 'tacint'),4), #additional parameters matrix
  MANIFESTVAR=matrix('merror',1), #measurement error
  T0VAR=diag(1e-3,4), #initial (t0) state uncertainty matrix (cholesky factor cov)
  T0MEANS=matrix(c('m1','m2','m3','taT0'),4), #initial state parameters
  MANIFESTMEANS=diag(0,1)) #measurement intercepts

#specify custom transforms
cm$pars$meanscale[cm$pars$param %in% 'p'] <- 100 #changed because this parameter had a large value
cm$pars$transform[cm$pars$param %in% 'tadynerrorsd'] <- 'log(1+exp(param))' #positive but not exponential transform
cm$pars$offset[cm$pars$param %in% 'm1'] <- ldat$X_obs.value[1] #we can estimate this starting value very well from data
cm$pars$transform[cm$pars$param %in% 'tadyn'] <- '-log(1+exp(param))' #negative but not exponential transform

cm$pars$indvarying<- FALSE #no variation over any parameters because single subject

#plot priors
plot(cm)

#specify custom state transition function (with reference to state and gradient vectors)
cm$gradient <- paste0( 
  '{
    vector[4] gradupd;
    real a = log(1+exp(state[4]));
    gradupd[1] = a * (state[2] - state[1]);
    gradupd[2] = state[1] * (PARS[2,1] - state[3]) - state[2];
    gradupd[3] = state[1] * state[2] - PARS[3,1] * state[3];
    gradupd[4] = state[4] * PARS[1,1] + PARS[4,1];
    gradient = gradupd;
  }')

#fit using ctsem
cftv <- ctStanFit(datalong = ldat,ctstanmodel = cm,
  optimcontrol=list(deoptim = F,isloops=2,isloopsize=500,issamples=500),optimize=T,verbose=1, nopriors = T, #remove this line for HMC instead of optimize + importance sample
  ukfspread=1e-1, 
  iter=300,chains=3)

summary(cftv)
par(mfrow=c(3,3))
ctStanPostPredict(cftv,diffsize = c(1,5),wait=FALSE)

#to retry with init using differential evolution
cfo <- optimstan(standata = cftv$standata,sm = cftv$stanmodel,deoptim = TRUE, decontrol = list(NP=40, steptol=50),verbose=1)
cftv$stanfit <- cfo




#sde specification using ctsem
cm <- ctModel(n.latent=3, n.TDpred = 0,n.manifest=1,type='stanct',
  manifestNames = 'X_obs.value',
  DRIFT=diag(-1,3),
  # DIFFUSION=diag(1e-5,3), #default is free parameter matrix
  LAMBDA=matrix(c(1,0,0),1),
  CINT=matrix(0,3,1),
  PARS = matrix(c('a','p','B'),3),
  MANIFESTVAR=matrix('merror',1),
  T0VAR=diag(1e-3,3), #this initial state cholesky cov needs to be fixed to 'something' when single subject
  T0MEANS=matrix(c('m1','m2','m3'),3),
  MANIFESTMEANS=diag(0,1))

cm$pars$transform[cm$pars$param %in% 'a'] <- 'exp(param)'
cm$pars$multiplier[cm$pars$matrix %in% 'DIFFUSION' & cm$pars$row == cm$pars$col] <- .1
cm$pars$meanscale[cm$pars$param %in% 'p'] <- 100
cm$pars$offset[cm$pars$param %in% 'm1'] <- ldat$X_obs.value[1]

cm$gradient <- paste0(
  '{
    vector[3] gradupd;
    gradupd[1] = PARS[1,1] * (state[2] - state[1]);
    gradupd[2] = state[1] * (PARS[2,1] - state[3]) - state[2];
    gradupd[3] = state[1] * state[2] - PARS[3,1] * state[3];
    gradient = gradupd;
  }')
    

#fit using ctsem
cf <- ctStanFit(datalong = ldat,ctstanmodel = cm,fit=T,
  optimize=T,verbose=1, nopriors = T,optimcontrol=list(deoptim=FALSE,isloops=2), #remove this line for sampling
  iter=300,chains=3)

summary(cf)
ctStanPostPredict(cf,diffsize = c(1,5),wait=FALSE)





#ode specification using ctsem
cm <- ctModel(n.latent=3, n.TDpred = 0,n.manifest=1,type='stanct',
  manifestNames = 'X_obs.value',
  DRIFT=diag(-1,3),
  DIFFUSION=diag(1e-8,3),
  LAMBDA=matrix(c(1,0,0),1),
  CINT=matrix(0,3,1),
  PARS = matrix(c('a','p','B'),3),
  MANIFESTVAR=matrix(.2,1),
  T0VAR=diag(1e-5,3),
  T0MEANS=matrix(c('m1','m2','m3'),3),
  MANIFESTMEANS=diag(0,1))

cm$pars$transform[cm$pars$param %in% 'a'] <- 'exp(param)'
cm$pars$meanscale[cm$pars$param %in% 'p'] <- 100
cm$pars$offset[cm$pars$param %in% 'm1'] <- ldat$X_obs.value[1]

cm$gradient <- paste0(
  '{
    vector[3] gradupd;
    gradupd[1] = PARS[1,1] * (state[2] - state[1]);
    gradupd[2] = state[1] * (PARS[2,1] - state[3]) - state[2];
    gradupd[3] = state[1] * state[2] - PARS[3,1] * state[3];
    gradient = gradupd;
  }')
    

#fit using ctsem
cf <- ctStanFit(datalong = ldat,ctstanmodel = cm,fit=T,optimize=T,verbose=1,nopriors = T,
  optimcontrol=list(deoptim=FALSE,isloops=2), #remove this line for sampling
  iter=300,chains=3)

summary(cf)
ctStanPostPredict(cf,diffsize = c(1,5),wait=F)

#try harder at optimizing
cfo <- optimstan(standata = cf$standata,cf$stanmodel,cores=3,verbose=1,issamples=500, isloops = 2, nopriors=TRUE,
  deoptim=T,decontrol=list(NP=40,steptol=30,c=0.05,CR=.5,strategy=6,p=.2))

cf$stanfit <- cfo

#update fit object to allow (simple) model changes without recompiling
cf <- ctStanUpdModel(cf,datalong = ldat,ctstanmodel = cm)