particleFilterLL_piecewise {pttstability}

run particle filter across piecewise data

Description

Calculates likelihoods across several segments of data - e.g. multiple plots from a single experiment. Requires several implicitely defined variables to run:

Usage

particleFilterLL_piecewise(
  param,
  N,
  y,
  libuse_y,
  smap_coefs,
  Euse,
  tuse,
  colpar = c(logit(1e-06), log(0.1)),
  nsmp = 1,
  lowerbound = -999,
  maxNuse = 512000
)

Arguments

param

parameters to be passed to parseparam0 function

N

number of particles

y

the time series to be analyzed

libuse_y

a matrix with two columns, specifying the start end end positions of segments within vector y

smap_coefs

a matrix of s-mapping coefficients

Euse

embedding dimension for the s-mapping analysis

tuse

theta for s-mapping analysis

colpar

parameters to be passed to the colfun0 - defaults to c(logit(1e-6), log(0.1))

nsmp

number of sample particle trajectories to return - defaults to 1

lowerbound

minimum accepted likelihood - used to automatically select number of particles. Defaults to -999

maxNuse

maximum number of particles to simulate - defaults to 512000

Value

results from particle filter - including mean estimates (Nest) and standard deviations (Nsd), across particles, and sample particle trajectories with (Nsmp) and without (Nsmp_noproc) process noise

Examples

# load data
data(dat)
# sort by index
dat = dat[order(dat$treatment, dat$number, dat$time),]


# make list of starting and ending positions for each replicate in the dat list
libmat<-NULL
trtmat<-data.frame(trt=as.character(sort(unique(dat$treatment))))
datnum<-1:nrow(dat)

for(i in 1:nrow(trtmat)) {
  ps1<-which(dat$treatment==trtmat$trt[i])
  replst<-sort(unique(dat$number[ps1]))

  for(j in 1:length(replst)) {
    ps2<-which(dat$number[ps1]==replst[j])
    libmat<-rbind(libmat, data.frame(trt=trtmat$trt[i], rep=replst[j],
      start=min(datnum[ps1][ps2]), end=max(datnum[ps1][ps2])))
  }
}

## run particle filter
# select treatment to analyse: enter either "LSA" or "LSP"
trtuse<-"HSP"
# extract library positions for treatment
libuse<-as.matrix(libmat[libmat$trt==trtuse,c("start", "end")])
# save abundance data to variable y
yps<-which(dat$treatment==trtuse)
y<-dat[,"Chlamydomonas.terricola"][yps]
libuse_y<-libuse-min(libuse)+1 # translate positions in dat to positions in y vector
y<-y/sd(y) # standardize to mean of one
timesteps<-dat$time[yps]

# get S-mapping parameters
sout<-data.frame(E = 2:4, theta = NA, RMSE = NA)
for(i in 1:nrow(sout)) {
  optout = optimize(f = function(x) {S_map_Sugihara1994(Y = y, E = sout$E[i],
      theta = x, lib = libuse_y)$RMSE}, interval = c(0,10))
  sout$theta[i] = optout$minimum
  sout$RMSE[i] = optout$objective
}
tuse<-sout$theta[which.min(sout$RMSE)] # find theta (nonlinerity) parameter
euse<-sout$E[which.min(sout$RMSE)] # find embedding dimension
spred<-S_map_Sugihara1994(Y = y, E = euse, theta = tuse, lib = libuse_y)

# set priors (log-transformed Beta_obs, Beta_proc1, and Beta_proc2)
minvUSE_edm<-c(log(0.001), log(0.001))  # lower limits
maxvUSE_edm<-c(log(2), log(2)) # upper limits


## Not run: 
  ## Run filter
  # Commented-out code: Install BayesianTools package if needed
  #install.packages("BayesianTools")
  set.seed(2343)
  require(BayesianTools)
  density_fun_USE_edm<-function(param) density_fun0(param = param,
    minv = minvUSE_edm, maxv=maxvUSE_edm)
  sampler_fun_USE_edm<-function(x) sampler_fun0(n = 1, minv = minvUSE_edm, maxv=maxvUSE_edm)
  prior_edm <- createPrior(density = density_fun_USE_edm, sampler = sampler_fun_USE_edm,
                         lower = minvUSE_edm, upper = maxvUSE_edm)
  niter<-5e3 # number of steps for the MCMC sampler
  N<-2e3 # number of particles
  smap_coefs<-process_scof(spred$C) # coefficients from s-mapping routine

  # likelihood and bayesian set-ups for EDM functions
  likelihood_EDM_piecewise_use<-function(x) {
    # default values for filter - see likelihood_EDM_piecewise documentation for details
    # note that colpar are set near zero because we expect no colonisation into a closed microcosm.
    likelihood_EDM_piecewise(param=x, y, libuse_y, smap_coefs, euse, tuse, N,
                             colpar = c(logit(1e-06), log(0.1)))
  }

  bayesianSetup_EDM <- createBayesianSetup(likelihood = likelihood_EDM_piecewise_use,
     prior = prior_edm)

  # run MCMC optimization (will take ~ 5 min)
  out_EDM <- runMCMC(bayesianSetup = bayesianSetup_EDM,
     settings = list(iterations=niter, consoleUpdates=20))
  burnin<-floor(niter/5) # burnin period
  plot(out_EDM, start=burnin) # plot MCMC chains
  gelmanDiagnostics(out_EDM, start=burnin) # calculate Gelman statistic
  summary(out_EDM, start=burnin) # coefficient summary

  ## extract abundance estimate from particle filter
  # use final estimate from MCMC chain
  smp_EDM<-(getSample(out_EDM, start=floor(niter/5)))
  tmp<-particleFilterLL_piecewise(param = smp_EDM[nrow(smp_EDM),], N=N, y = y, libuse_y = libuse_y,
                                  smap_coefs = smap_coefs, Euse = euse, tuse = tuse)
  # mean estimated abundance
  simout<-tmp$Nest
  # sd estimated abundance
  sdout<-tmp$Nsd
  # sample from true particle trajectory
  simout_smp<-tmp$Nsmp
  # sample from true particle trajectory pre-process noise
  simout_smp_noproc<-tmp$Nsmp_noproc

  plot(timesteps, simout, xlab="Time", ylab="Abundance")
  abline(h=0, lty=3)

## End(Not run)