The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
bayesSSM is an R package offering a set of tools for performing
Bayesian inference in state-space models (SSMs). It implements the
Particle Marginal Metropolis-Hastings (PMMH) in the main function
pmmh
for Bayesian inference in SSMs.
While there are several alternative packages available for performing Particle MCMC bayesSSM is designed to be simple and easy to use. It was developed as a procrastination task during my Master’s thesis about Particle MCMC, since I was implementing everything from scratch anyway. Everything is written in R, so performance is not the best.
You can install the latest stable version of bayesSSM from CRAN with:
install.packages("bayesSSM")or the development version from GitHub with:
# install.packages("pak")
pak::pak("BjarkeHautop/bayesSSM")Consider the following SSM:
\[ \begin{aligned} X_1 &\sim N(0,1) \\ X_t&=\phi X_{t-1}+\sin(X_{t-1})+\sigma_x V_t, \quad V_t \sim N(0,1) \\ Y_t&=X_t+\sigma_y W_t, \quad W_t \sim N(0, \, 1). \end{aligned} \]
Let’s first simulate some data from this model with \(\phi = 0.8\), \(\sigma_x = 1\), and \(\sigma_y = 0.5\).
t_val <- 20
phi_val <- 0.8
sigma_x_val <- 1
sigma_y_val <- 0.5
x <- numeric(t_val)
y <- numeric(t_val)
x[1] <- rnorm(1, mean = 0, sd = sigma_x_val)
y[1] <- rnorm(1, mean = x[1], sd = sigma_y_val)
for (t in 2:t_val) {
x[t] <- phi_val * x[t - 1] + sin(x[t - 1]) + rnorm(1,
mean = 0,
sd = sigma_x_val
)
y[t] <- x[t] + rnorm(1, mean = 0, sd = sigma_y_val)
}We define the priors for our model as follows:
\[ \begin{aligned} \phi &\sim \text{Uniform}(0,1), \\ \sigma_x &\sim \text{Exp}(1), \\ \sigma_y &\sim \text{Exp}(1). \end{aligned} \]
We can use pmmh to perform Bayesian inference on this
model. To use pmmh we need to define the functions for the
SSM and the priors. The functions init_fn,
transition_fn should be functions that simulates the latent
states. They must contain the argument particles, which is
a vector of particles, and can contain any other arguments. The function
log_likelihood_fn should be a function that calculates the
log-likelihood of the observed data given the latent state variables. It
must contain the arguments y and
particles.
The priors for the parameters must be defined as log-prior functions.
Every parameter from init_fn, transition_fn,
and log_likelihood_fn must have a corresponding log-prior
function.
init_fn <- function(particles) {
stats::rnorm(particles, mean = 0, sd = 1)
}
transition_fn <- function(particles, phi, sigma_x) {
phi * particles + sin(particles) +
stats::rnorm(length(particles), mean = 0, sd = sigma_x)
}
log_likelihood_fn <- function(y, particles, sigma_y) {
stats::dnorm(y, mean = particles, sd = sigma_y, log = TRUE)
}
log_prior_phi <- function(phi) {
stats::dunif(phi, min = 0, max = 1, log = TRUE)
}
log_prior_sigma_x <- function(sigma) {
stats::dexp(sigma, rate = 1, log = TRUE)
}
log_prior_sigma_y <- function(sigma) {
stats::dexp(sigma, rate = 1, log = TRUE)
}
log_priors <- list(
phi = log_prior_phi,
sigma_x = log_prior_sigma_x,
sigma_y = log_prior_sigma_y
)Now we can run the PMMH algorithm using the pmmh
function. We run 2 chains for 200 MCMC samples with a burn-in of 10. We
also modify the tuning to only use 200 pilot samples and a burn-in of
10. In practice you would want to run it for a much larger number of
samples.
library(bayesSSM)
result <- pmmh(
y = y,
m = 500, # number of MCMC samples
init_fn = init_fn,
transition_fn = transition_fn,
log_likelihood_fn = log_likelihood_fn,
log_priors = log_priors,
pilot_init_params = list(
c(phi = 0.4, sigma_x = 0.4, sigma_y = 0.4),
c(phi = 0.8, sigma_x = 0.8, sigma_y = 0.8)
),
burn_in = 50,
num_chains = 2,
seed = 1405,
tune_control = default_tune_control(pilot_m = 200, pilot_burn_in = 10)
)
#> Running chain 1...
#> Running pilot chain for tuning...
#> Using 50 particles for PMMH:
#> Running particle MCMC chain with tuned settings...
#> Running chain 2...
#> Running pilot chain for tuning...
#> Using 53 particles for PMMH:
#> Running particle MCMC chain with tuned settings...
#> PMMH Results Summary:
#> Parameter Mean SD Median CI Lower.2.5% CI Upper.97.5% ESS Rhat
#> phi 0.66 0.15 0.67 0.38 0.94 25 1.036
#> sigma_x 0.84 0.33 0.88 0.17 1.43 11 1.169
#> sigma_y 0.63 0.38 0.56 0.07 1.51 16 1.170
#> Warning in pmmh(y = y, m = 500, init_fn = init_fn, transition_fn =
#> transition_fn, : Some ESS values are below 400, indicating poor mixing.
#> Consider running the chains for more iterations.
#> Warning in pmmh(y = y, m = 500, init_fn = init_fn, transition_fn = transition_fn, :
#> Some Rhat values are above 1.01, indicating that the chains have not converged.
#> Consider running the chains for more iterations and/or increase burn_in.We get convergence warnings as expected due to the small number of iterations.
A state-space model (SSM) has the structure given in the following directed acyclic graph (DAG):
The core function, pmmh, implements the Particle
Marginal Metropolis-Hastings, which is an algorithm that first generates
a set of \(N\) particles to approximate
the likelihood and then uses this approximation in the acceptance
probability. The implementation automatically tunes the number of
particles and the proposal distribution for the parameters.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.