We will show how to fit the following SSM model using the
bayesSSM package:
X0∼N(0,1)Xt=ϕXt−1+sin(Xt−1)+σxVt,Vt∼N(0,1),t≥1Yt=Xt+σyWt,Wt∼N(0,1),t≥1 that is Xt is a latent state and Yt is an observed value. The parameters of the model are ϕ, σx, and σy.
First, we will simulate some data from this model:
set.seed(1405)
t_val <- 20
phi <- 0.8
sigma_x <- 1
sigma_y <- 0.5
init_state <- rnorm(1, mean = 0, sd = 1)
x <- numeric(t_val)
y <- numeric(t_val)
x[1] <- phi * init_state + sin(init_state) +
rnorm(1, mean = 0, sd = sigma_x)
y[1] <- x[1] + rnorm(1, mean = 0, sd = sigma_y)
for (t in 2:t_val) {
x[t] <- phi * x[t - 1] + sin(x[t - 1]) + rnorm(1, mean = 0, sd = sigma_x)
y[t] <- x[t] + rnorm(1, mean = 0, sd = sigma_y)
}
x <- c(init_state, x)Let’s visualize the data:
ggplot() +
geom_line(aes(x = 0:t_val, y = x), color = "blue", linewidth = 1) + # Latent
geom_point(aes(x = 1:t_val, y = y), color = "red", size = 2) + # Observed
labs(
title = "Simulated Data: Latent State and Observations",
x = "Time",
y = "Value",
caption = "Blue line: Latent state (x), Red points: Observed values (y)"
) +
theme_minimal()
To fit the model using pmmh we need to specify the
initialization, transition, and log-likelihood functions. The
initialization function init_fn must take the argument
num_particles and return a vector or matrix of particles.
The transition function transition_fn must take the
argument particles and return a vector or matrix of
particles. The log-likelihood function log_likelihood_fn
must take the arguments y, particles, and
return a vector of log-likelihood values.
All the functions can take model-specific parameters as arguments.
Time-dependency can be implemented by giving a t argument
in transition_fn and log_likelihood_fn.
init_fn <- function(num_particles) {
rnorm(num_particles, mean = 0, sd = 1)
}
transition_fn <- function(particles, phi, sigma_x) {
phi * particles + sin(particles) +
rnorm(length(particles), mean = 0, sd = sigma_x)
}
log_likelihood_fn <- function(y, particles, sigma_y) {
dnorm(y, mean = particles, sd = sigma_y, log = TRUE)
}Since we are interested in Bayesian inference, we need to specify the
priors for our parameters. We will use a normal prior for ϕ and exponential priors for σx and σy. pmmh needs the
priors to be specified on the log-scale and takes the priors as a list
of functions.
log_prior_phi <- function(phi) {
dunif(phi, min = 0, max = 1, log = TRUE)
}
log_prior_sigma_x <- function(sigma) {
dexp(sigma, rate = 1, log = TRUE)
}
log_prior_sigma_y <- function(sigma) {
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
)The pmmh function automatically tunes the number of
particles and proposal distribution for the parameters. The tuning can
be modified by the the function default_tune_control.
We fit 2 chains with m=1000 iterations for each, with a burn_in of 500. We also modify the tuning to only use a pilot run of 100 iterations and 10 burn-in iterations. In practice you should run more iterations and chains. To improve sampling we specify that proposals for σx and σy should be on the log-scale.
result <- pmmh(
y = y,
m = 1000,
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 = 500,
num_chains = 2,
seed = 1405,
param_transform = list(
phi = "logit",
sigma_x = "log",
sigma_y = "log"
),
tune_control = default_tune_control(pilot_m = 100, pilot_burn_in = 10),
verbose = TRUE
)
#> Running chain 1...
#> Running pilot chain for tuning...
#> Pilot chain posterior mean:
#> phi sigma_x sigma_y
#> 0.8504781 0.8465155 0.9828058
#> Pilot chain posterior covariance (on transformed space):
#> phi sigma_x sigma_y
#> phi 0.002853161 -0.001659855 0.005267221
#> sigma_x -0.001659855 0.030056091 -0.019598925
#> sigma_y 0.005267221 -0.019598925 0.025940662
#> Using 50 particles for PMMH:
#> Running Particle MCMC chain with tuned settings...
#> Running chain 2...
#> Running pilot chain for tuning...
#> Pilot chain posterior mean:
#> phi sigma_x sigma_y
#> 0.6443472 1.1777341 0.4654751
#> Pilot chain posterior covariance (on transformed space):
#> phi sigma_x sigma_y
#> phi 4.532081e-03 -0.011399852 -3.544264e-05
#> sigma_x -1.139985e-02 0.050376412 -1.354051e-03
#> sigma_y -3.544264e-05 -0.001354051 1.003767e-02
#> Using 343 particles for PMMH:
#> Running Particle MCMC chain with tuned settings...
#> PMMH Results Summary:
#> Parameter Mean SD Median 2.5% 97.5% ESS Rhat
#> phi 0.95 0.10 1.00 0.68 1.00 2 2.080
#> sigma_x 0.98 0.40 1.04 0.26 1.67 1 1.426
#> sigma_y 0.69 0.53 0.51 0.05 1.57 1 2.584
#> Warning in pmmh(y = y, m = 1000, 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 = 1000, 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 see that the chains gives convergence issues, indicating that we should run it for more iterations, but we ignore this issue in this Vignette.
It automatically prints data frame summarizing the results, which can
be printed from any pmmh_output object by calling
print.
print(result)
#> PMMH Results Summary:
#> Parameter Mean SD Median 2.5% 97.5% ESS Rhat
#> phi 0.95 0.10 1.00 0.68 1.00 2 2.080
#> sigma_x 0.98 0.40 1.04 0.26 1.67 1 1.426
#> sigma_y 0.69 0.53 0.51 0.05 1.57 1 2.584The chains are saved as theta_chain
Let’s collect the chains for phi from the chains and visualize the densities
ggplot(chains, aes(x = phi, fill = factor(chain))) +
geom_density(alpha = 0.5) +
labs(
title = "Density plot of phi chains",
x = "Value",
y = "Density",
fill = "Chain"
) +
theme_minimal()
We have now fitted a simple SSM model using bayesSSM.
Feel free to explore the package further and try out different
models.