The Particle Marginal Metropolis Hastings (PMMH) algorithm can be used to learn parameters of a Dynamic Generalised Linear Model (DGLM).

The algorithm uses an unbiased estimate of the marginal likelihood, \(p(y_{1:T} | x_{0:T}, V, W)\) in a Metropolis-Hastings Algorithm. The pseudo marginal likelihood can be calculated using the Particle Filter.

Firstly define a DGLM:

scala> import com.github.jonnylaw.dlm._
<console>:12: warning: Unused import
       import com.github.jonnylaw.dlm._
                                      ^
import com.github.jonnylaw.dlm._

scala> val model = Dlm.polynomial(1)
model: com.github.jonnylaw.dlm.Dlm = Dlm(com.github.jonnylaw.dlm.Dlm$$$Lambda$6788/1540206298@3523d50b,com.github.jonnylaw.dlm.Dlm$$$Lambda$6789/1989125347@26334918)

scala> val dglm = Dglm.poisson(model)
dglm: com.github.jonnylaw.dlm.Dglm = Dglm(com.github.jonnylaw.dlm.Dglm$$$Lambda$7171/2056819218@77b45cca,com.github.jonnylaw.dlm.Dlm$$$Lambda$6788/1540206298@3523d50b,com.github.jonnylaw.dlm.Dlm$$$Lambda$6789/1989125347@26334918,com.github.jonnylaw.dlm.Dglm$$$Lambda$7172/1171156020@4982582c,com.github.jonnylaw.dlm.Dglm$$$Lambda$7173/490903945@7f9b32e4)

Assume we have some observations of a process which can be modelled using the specified model:

val observations = // some observations

Then we need to define the proposal function, def proposal: Parameters => Rand[Parameters]. A symmetric proposal distribution is implemented in the Metropolis object. This proposal distribution is simply a Normal centered at the previous value of the parameters. In the case of the non-negative parameters, \(V\), \(W\) and \(C_0\) the parameter is multiplied by the exponential of a Normal random number.

The function accepts a parameter delta, a single Double value which controls the variance of the Gaussian proposal distribution:

scala> def proposal = Metropolis.symmetricProposal(0.05) _
proposal: com.github.jonnylaw.dlm.DlmParameters => breeze.stats.distributions.Rand[com.github.jonnylaw.dlm.DlmParameters]

Then we need to specify a prior distribution for the parameters, which encapsulates our prior beliefs about the parameters.

def prior(p: Dlm.Parameters): Double = ???

Then we must specify the initial value of the parameters to start the MCMC, this can be a sample from the prior distribution. The PMMH algorithm is available to use in the Bayesian DLM Scala library, using the bootstrap particle filter to estimate the likelihood:

val iters = MetropolisHastings.pmmh(
  dglm, 
  observations,
  proposal,
  prior,
  initP
  n = 500)

With this specification, there are two tuning parameters in the PMMH algorithm, the variance of the Gaussian proposal distribution and the number of particles in the Particle Filter. As the number of particles increases so does the accuracy of the likelihood estimate, however this also increases computational time. The variance of the proposal distribution should be such that the proportion of accepted moves in the pmmh algorithm is approximately one third. This will guarantee that the MCMC has explored the full parameter posterior.