Fitting Markov Modulated Point Process Models

Owen G. Ward.

library(ppdiag)
library(rstan)
library(cmdstanr)
library(tidyverse)
library(bayesplot)

Functions for fitting standard Hawkes and Poisson point processes to data are included in ppdiag. However, currently, we do not include fitting algorithms for Markov modulated point processes, as these rely on the use of rstan or (more recently), cmdstanr to perform Bayesian inference for the model parameters.

Here we provide instructions on how to fit these models so that they can easily be used in conjunction with the diagnostic tools of ppdiag.

Installing RStan or Cmdstanr

• Instructions for installing RStan can be found at https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started
• Instructions for installing CmdStanR can be found at https://mc-stan.org/cmdstanr/articles/cmdstanr.html. As mentioned in this link there are some differences between RStan and CmdStanR. In our experience, CmdStanR does seem to be maybe more reliable and less prone to crashing.
• If you are having trouble installing or using Stan there is lots of great help available at https://discourse.mc-stan.org/, where someone else has probably already answered your question!

The Stan code for MMPP/MMHP

We first include the stan code to fit each of these Markov modulated models to a temporal point process.

mmpp_stan_code <- [2300 chars quoted with '"']

mmhp_stan_code <- [6866 chars quoted with '"']

Fitting these models in RStan

To demonstrate how to use these models, we will simulate some data from each of these models and fit the included stan models.

Q <- matrix(c(-0.4, 0.4, 0.2, -0.2), ncol = 2, byrow = TRUE)
mmpp_obj <- pp_mmpp(Q = Q, lambda0 = 1, c = 1.5, delta = c(1/3, 2/3))

sim_mmpp <- pp_simulate(mmpp_obj, n = 50)

We also show how to run these models using rstan for the MMPP example.

mmpp_data <- list(num_events = length(sim_mmpp$events) - 1, 
                  # as first event is start time (0)
                  time_matrix = diff(c(sim_mmpp$events, sim_mmpp$end))
                  #interevent arrival time
                  )

mmpp_rstan <- stan(model_code = mmpp_stan_code, 
                   data = mmpp_data,
                   chains = 2)
mmpp_sim <- rstan::extract(mmpp_rstan)

The fit of this model can be evaluated in several ways, including some of the tools in the bayesplot package, which are important to ensure the fit is reasonable.

bayesplot::mcmc_hist(as.matrix(mmpp_rstan), pars = c("lambda0", "c"))

To use the results of this fit with the functions of ppdiag, we want to extract an estimate of each of the parameters. This can be achieved using the posterior mean.

The desired posterior quantities can then be extracted from this object.

mmpp_post <- lapply(mmpp_sim, mean)

Similarly, we can fit the MMHP model to simulated data also.

Q <- matrix(c(-0.4, 0.4, 0.2, -0.2), ncol = 2, byrow = TRUE)
mmhp_obj <- pp_mmhp(Q = Q, lambda0 = 0.5, lambda1 = 1.5,
                    alpha = 0.5, beta = 0.75, delta = c(1/3, 2/3))

sim_mmhp <- pp_simulate(mmhp_obj, n = 50)


mmhp_data <- list(num_events = length(sim_mmhp$events) - 1, 
                  # as first event is start time (0)
                  time_matrix = diff(c(sim_mmhp$events, sim_mmhp$end))
                  #interevent arrival time
                  )

mmhp_rstan <- stan(model_code = mmhp_stan_code, 
                   data = mmhp_data,
                   chains = 2)
mmhp_sim <- rstan::extract(mmhp_rstan)
mmhp_post <- lapply(mmhp_sim, mean)

bayesplot::mcmc_hist(as.matrix(mmhp_rstan),
                     pars = c("lambda0", "lambda1", "alpha", "beta"))

Fitting these models in Cmdstanr

We can also fit these models using cmdstanr. We include the code to fit these examples below (which is not run here).

mmpp_file <- write_stan_file(mmpp_stan_code)
mmpp_stan <- cmdstan_model(stan_file = mmpp_file)

fit_mmpp <- mmpp_stan$sample(data = mmpp_data,
                             seed = 123,
                             chains = 4,
                             parallel_chains = 4,
                             refresh = 500)

Similarly, we can fit the MMHP model in cmdstanr.

mmhp_file <- write_stan_file(mmhp_stan_code)
mmhp_stan <- cmdstan_model(stan_file = mmhp_file)

fit_mmhp <- mmhp_stan$sample(data = mmhp_data,
                             seed = 123,
                             chains = 4,
                             parallel_chains = 4,
                             refresh = 500)

Using the Stan output with ppdiag

Having obtained estimates from fitting these stan models, we can now pass these parameter estimates into ppdiag to evaluate model fit.

mmpp_fit_obj <- pp_mmpp(lambda0 = mmpp_post$lambda0,
                       c = mmpp_post$c,
                       Q = matrix( c(-mmpp_post$q1,
                                     mmpp_post$q1,
                                     mmpp_post$q2,
                                     -mmpp_post$q2),
                                   nrow = 2, ncol = 2, byrow = T) )


pp_diag(mmpp_fit_obj, events = sim_mmpp$events)

mmhp_fit_obj <- pp_mmhp(lambda0 = mmhp_post$lambda0,
                        lambda1 = mmhp_post$lambda1,
                        alpha = mmhp_post$alpha,
                        beta = mmhp_post$beta,
                        Q = matrix( c(-mmhp_post$q1,
                                      mmhp_post$q1,
                                      mmhp_post$q2,
                                      -mmhp_post$q2),
                                 nrow = 2, ncol = 2, byrow = T) )

pp_diag(mmhp_fit_obj, events = sim_mmhp$events)