SMC2_ABC: Generate posterior samples of static and state parameters.
In AnthonyEbert/StateSpaceInference:

Description Usage Arguments Details Examples

View source: R/SMC2_ABC.R

Generate posterior samples of static and state parameters.

SMC2_ABC(
  prior_sample,
  dprior,
  loss,
  loss_args,
  Ntheta,
  Nx,
  pacc,
  dtp = 1,
  ESS_threshold = 0.1,
  eps = NULL,
  cl = NULL,
  TT,
  trans = function(x, trans_args) {     I(x) },
  invtrans = function(x, trans_args) {     I(x) },
  resample_times = NA,
  trans_args = list(),
  cov_coef = 1,
  acceptance_correction = function(x) {     1 }
)

`prior_sample`	Matrix of prior samples of static parameters. The number of rows is the number of samples, the number of components of the parameter is equal to the number of columns.
`dprior`	Function which evaluates the prior density. This function should have one argument taking the entire vector of static parameters.
`loss`	Function which takes as input values of static and state parameters and returns a distance and an updated value for the state parameter, see Details.
`loss_args`	List argument which is passed to the `loss` function. See Details.
`Ntheta`	Positive number representing the number of static parameter proposals.
`Nx`	Positive number representing the number of state parameter proposals for each static parameter proposal.
`pacc`	Number between 0 and 1. The threshold for calculating the threshold for accepting state parameter proposals.
`dtp`	Positive number. Time between time steps.
`ESS_threshold`	Positive number. Effective sample size (ESS) threshold, when ESS drops below this number multiplied by `Ntheta` then static parameter proposals are replenished.
`eps`	Positive valued numeric vector of length `TT`, optional. This vector is usually generated automatically from `pacc`, but if this argument is used then `pacc` is ignored.
`cl`	An object of class "cluster" to use `parLapply` internally or the string "mclapply" to use `mclapply` internally. A NULL input (the default) will use `lapply` internally, a single core.
`TT`	Positive integer. Number of time steps.
`trans`	Function to transform parameters for resampling purposes, see Details.
`invtrans`	Inverse function of `trans` to invert the transformation.
`resample_times`	Numeric vector of replenishment times, optional.
`trans_args`	List argument which is passed to the `trans` and `invtrans` function.
`cov_coef`	Positive number to
`acceptance_correction`	Positive number. Correction to acceptance probability for `trans` function.

The loss argument is a function with the following arguments: x, theta, time1, time2, inp. The argument x is the current state parameter proposal, theta is the static parameter proposal, time1 and time2 are the start and end times for the update step, and inp takes a list of any other inputs to the function. This function should return list of length two with a distance and an updated value for the state parameter. See the example.

## Not run: 
library(parallel)
library(StateSpaceInference)
library(ggplot2)
library(ggalt)

cl <- makeCluster(parallel::detectCores() - 1)
#cl <- "mclapply"
#cl <- NULL

TT <- 10
true_theta <- c(0.25, 0.5)
lower <- 0
upper <- 3.5
sd_t <- 1
init <- min(rgamma(1, 100, 100), upper - 1)
a_logit <- 0.9
dist_coef <- 0.5
true_states <-cumsum(rnorm(TT))

lambda_fun <- stepfun(seq(1, TT - 1, by = 1), y = true_states)

generator <- function(TT, true_states, theta){
y <- list()
for(i in 1:TT){
  y[[i]] <- as.numeric(sn::rsn(1e3, dp = sn::cp2dp(cp = c(true_states[i], theta), "SN")))
}
return(y)
}

y <- generator(TT, true_states, true_theta)


inp <- list(
  lower = lower,
  upper = upper,
  sd_t = sd_t,
  a_logit = a_logit,
  y = y
)

loss <- function(x, theta, time1, time2, inp){
x <- as.numeric(x[1])

if(gtools::invalid(x)){
  x <- rnorm(1)
} else {
  x <- rnorm(1, mean = x)
}

y <- sn::rsn(1e1, dp = sn::cp2dp(cp = c(x, theta), "SN"))

ss_obs <- c(mean(inp$y[[time2]]), sd(inp$y[[time2]]), e1071::skewness(inp$y[[time2]]))
ss_sim <- c(mean(y), sd(y), e1071::skewness(y))


return(list(distance = sqrt(sum((ss_obs - ss_sim)^2)), x = x))
}


Ntheta = 200
Nx = 10
pacc = 0.5

lower_theta <- c(0.1, 0.2)
upper_theta <- c(0.5, 0.8)

prior_sample <- data.frame(theta1 = runif(Ntheta, lower_theta[1], upper_theta[1]), theta2 = runif(Ntheta, lower_theta[2], upper_theta[2]))

prior_sample <- as.matrix(prior_sample, ncol = 2)

trans <- function(x, trans_args){
  theta1 <- log(x[,1])
  theta2 <- log(x[,2])
  return(cbind(theta1, theta2))
}

invtrans <- function(x, trans_args){
  theta1 <- exp(x[,1])
  theta2 <- exp(x[,2])
  return(cbind(theta1, theta2))
}

full_list <- SMC2_ABC(prior_sample, dprior = function(x){dunif(x[1], 0.1, 0.5)*dunif(x[2],0.2,0.8)}, loss, loss_args = inp, Ntheta = Ntheta, Nx = Nx, pacc = pacc, cl = cl, dt = 1, ESS_threshold = 0.5, TT = TT, trans = trans, invtrans = invtrans, cov_coef = 0.25^2)

state_df <- get_state(full_list)

state_df$state <- true_states

theta_df <- get_parameter(full_list)


ggplot(state_df) + aes(x = time, y = state, ymin = lower, ymax = upper) + geom_step() + geom_ribbon(alpha = 0.2, stat = "stepribbon", fill = "red") + geom_step(mapping = aes(x = time, y = med), col = "red") + ggthemes::theme_base() + scale_y_continuous(expand = c(0, 0)) + scale_x_continuous(expand = c(0, 0))

ggplot(theta_df[which(theta_df$time %% 5 == 0),]) + aes(x = value, weights = weight, col = factor(time)) + geom_density() + facet_wrap(~parameter, scales = "free")

## End(Not run)