4
$\begingroup$

I am following the Hidden Markov Models guide text for Time Series An Introduction Using R (Walter Zucchini). Chapter 7. Bayesian inference for Poisson-hidden Markov models. Specifically in section 7.2 talk about how to select the right model from the Bayes factor between two models.

Then use a code in rstan to estimate the Poisson-Hidden Markov model parameters. For the base data Earthquakes: Number of major earthquakes (magnitude 7 or greater) in the World, 1900-2006.

library(rstan)

dat <- read.table("http://www.hmms-for-time-series.de/second/data/earthquakes.txt")
stan.data <- list(n=dim(dat)[1], m=2, x=dat$V2)

PHMM <- '
data {
  int<lower=0> n; // length of the time series
  int<lower=0> x[n]; // data
  int<lower=1> m; // number of states 
}

parameters{
  simplex[m] Gamma[m]; // tpm
  positive_ordered[m] lambda; // mean of poisson - ordered
}  

model{
  vector[m] log_Gamma_tr[m]; // log, transposed tpm 
  vector[m] lp; // for forward variables
  vector[m] lp_p1; // for forward variables

  lambda ~ gamma(0.1, 0.01); // assigning exchangeable priors 
  //(lambdas´s are ordered for sampling purposes)

  // transposing tpm and taking the log of each entry
  for(i in 1:m)
    for(j in 1:m)
      log_Gamma_tr[j, i] = log(Gamma[i, j]);

  lp = rep_vector(-log(m), m); // 

    for(i in 1:n) {
      for(j in 1:m)
        lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp) + poisson_lpmf(x[i] | lambda[j]); 

      lp = lp_p1;
    }

  target += log_sum_exp(lp);
}'

The output is as follows:

modelo <- stan(model_code = PHMM, data = stan.data, iter = 1000, chains = 1)
print(modelo,digits_summary = 3)
Inference for Stan model: f4fa370704db17d55928f5a217f01f5d.
1 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=500.

               mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff  Rhat
Gamma[1,1]    0.905   0.003 0.049    0.786    0.881    0.912    0.939    0.973   280 1.004
Gamma[1,2]    0.095   0.003 0.049    0.027    0.061    0.088    0.119    0.214   280 1.004
Gamma[2,1]    0.142   0.003 0.068    0.032    0.093    0.135    0.178    0.308   501 0.998
Gamma[2,2]    0.858   0.003 0.068    0.692    0.822    0.865    0.907    0.968   501 0.998
lambda[1]    15.068   0.053 0.953   13.121   14.497   15.184   15.736   16.641   324 0.998
lambda[2]    25.611   0.082 1.450   22.948   24.675   25.657   26.504   28.830   310 0.998
lp__       -350.196   0.108 1.666 -354.707 -350.908 -349.797 -348.959 -348.187   237 0.998

Samples were drawn using NUTS(diag_e) at Thu Feb 07 14:24:17 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

7.2 Bayesian estimation of the number of states.

In the Bayesian approach to model selection, the number of states, $m$, is a parameter whose value is assessed from its posterior distribution, $p(m | x^{(T)})$. Computing this posterior distribution is, however, not an easy problem; indeed, it has been described as `notoriously difficult to calculate' (Scott, James and Sugar, 2005). Using $p$ as a general symbol for probability mass or density functions, one has

$$p(m | x^{(T)}) = p(m) p(x^{(T)} | m) / p(x^{(T)}) \propto p(m) p(x^{(T)} | m) \ \ \ \ \ (7.3)$$

where $p(x^{(T)} | m)$ is called the integrated likelihood. If only two models are being compared, the posterior odds are equal to the product of the `Bayes factor' and the prior odds:

$$\frac{p(m_2 | x^{(t)})}{p(m_1 | x^{(t)})} = \frac{p(x^{(T)} | m_2)}{p(x^{(T)} | m_1)} \times \frac{p(m_2)}{p(m_1)} \ \ \ \ \ (7.4)$$

7.2.1 Use of the integrated likelihood

In order to use (7.3) or (7.4) we need to estimate the integrated likelihood

$$p(x^{(T)} | m) = \int p(\theta_m, x^{(T)} | m) d \theta_m \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \int p(x^{(T)} | m, \theta_m) p(\theta_m | m) d \theta_m$$

One way of doing so would be to simulate from $p(\theta_m | m)$, the prior distribution of the parameters $\theta_m$ of the m-state model. But it is convenient and especially if the prior is diffuse and more efficient to use a method that requires instead a sample from the posterior distribution, $p(\theta_m | x(T), m)$. Such a method is as follows.

Write the integrated likelihood as

$$\int p(x^{(T)} | m, \theta_m) \frac{p(\theta_m | m)}{p^*(\theta_m)} p*(\theta_m) d \theta_m$$

that is, write it in a form suitable for the use of a sample from some convenient density $p^*(\theta_m)$ for the parameters $\theta_m$. Since we have available a sample $\theta_m^{(j)} \ (j=1,2,...,B)$ from the posterior distribution, we can use that sample; that is, we can take $p^*(\theta_m) = p(\theta_m | x^{(T)}, m)$. Newton and Raftery (1994) therefore suggest inter alia that the integrated likelihood can be estimated by

$$\hat{I} = \sum_{j=1}^B w_j \ p(x^{(T)} | m, \theta_m^{(j)}) \Big/ \sum_{j=1}^B w_j, \ \ \text{where} \ \ w_j = \frac{p(\theta_m^{(j)} | m)}{p(\theta_m^{(j)} | x^{(T)}, m)} \ \ \ \ \ \ (7.5)$$

After some manipulation this simplies to the harmonic mean of the likelihood values of a sample from the posterior,

$$\hat{I} = \left(B^{-1} \sum_{j=1}^B \left( p(x^{(T)} | m, \theta_m^{(j)}) \right)^{-1} \right); \ \ \ \ \ \ (7.6)$$

In section 7.2.2 Zucchini, indicates that another and simpler way to select the best model by parallel sampling.

The truth is that I have no idea how to elaborate the code in R, which allows me to calculate the factor of bays either by method 1. using the integrated likelihood or by method 2 using parallel sampling and thus choosing the most suitable model. Any suggestion is well-received.

New information

enter image description here

A hidden Markov model $\lbrace X_t : t \in \mathbb{N} \rbrace$ is a particular kind of dependent mixture. With $X^{(t)}$ and $C^{(t)}$ representing the histories from time $1$ to time $t$, one can summarize the simplest model of this kind by:

$$P_r (C_t | C^{(t-1)}) = P_r (C_t | C_{t-1}), \ \ \ t = 2,3,... \ \ \ \ \ (2.1)$$

$$P_r(X_t | X^{(t-1)}, C^{(t)}) = P_r(X_t | C_t) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2.2)$$

The model consists of two parts: firstly, an unobserved 'parameter process' $\lbrace C_t: t = 1,2,... \rbrace$ satisfying the Markov property; and secondly, the 'state-dependent process' $\brace X_t: t = 1,2,... \rbrace$, in which the distribution of $X_t$ depends only on the current state $C_t$ and not on previous states or observations. This structure is represented by the directed graph in Figure 2.2.

$\endgroup$
2
  • $\begingroup$ If there are some excerpts from the book, they should be clearly indicated as quotes. $\endgroup$
    – Xi'an
    Feb 11, 2019 at 9:06
  • $\begingroup$ And also please include the model as a mathematical formula, as readers cannot be expected to read Zucchini´s book or translate the rstan code. $\endgroup$
    – Xi'an
    Feb 11, 2019 at 9:42

2 Answers 2

3
+25
$\begingroup$

The method of Newton and Raftery (1994) is also called the harmonic mean estimator of the marginal/integrated likelihood and it has been dubbed the "worst Monte Carlo method ever" by Radford Neal (U Toronto) with very strong arguments that it can essentially produce results with no numerical connections to the truth, by arbitrary orders of magnitude. I have also written many warnings against using it and discussed methods that avoid the infinite variance danger.

I have no clear idea what parallel sampling stands for in the field of computational statistics. If W. Zucchini means parallel tempering this is a rather delicate technique that proves difficult to calibrate. If he means the method advocated by Peter Congdon, we demonstrated it is invalid, being based on a confusion between model-based posteriors and joint pseudo-posteriors.

$\endgroup$
0
$\begingroup$

A possible solution to this problem is to use bridging sampling, in the following link there is a complete tutorial of your application. In addition, I quote from this article the following conclusions: "The bridge sampling estimator is superior to the naive Monte Carlo estimator, theimportance sampling estimator, and the generalized harmonic mean estimator for severalreasons. First, Meng and Wong (1996) showed that, among the four estimators discussedin this article, the bridge sampler presented in this article minimizes the mean-squarederror because it uses the optimal bridge function. Second, in bridge sampling, choosing asuitable proposal distribution is much easier than choosing a suitable importance densityfor the importance sampling estimator or the generalized harmonic mean estimator becausebridge sampling is more robust to the tail behavior of the proposal distribution relativeto the posterior distribution." Gronau, et al. (2017) Assume two models, Poisson-HMM of two states and a Poisson-HMM of three states.

model_H0 <- stan(model_code = PHMM, 
                 data = list(n=dim(dat)[1], m=2, x=dat$V2), 
                 iter = 1000, chains = 1)

model_H1 <- stan(model_code = PHMM, 
               data = list(n=dim(dat)[1], m=3, x=dat$V2), 
               iter = 1000, chains = 1)

Use the bridgesampling library that has implemented the method to be used.

library(bridgesampling)
bridge_H0 <- bridge_sampler(model_H0)
bridge_H1 <- bridge_sampler(model_H1)
print(bridge_H0)
Bridge sampling estimate of the log marginal likelihood: -350.8634
Estimate obtained in 7 iteration(s) via method "normal".
print(bridge_H1)
Bridge sampling estimate of the log marginal likelihood: -348.5192
Estimate obtained in 6 iteration(s) via method "normal".

compute approximate percentage errors

error_measures(bridge_H0)$percentage
[1] "1.44%"
error_measures(bridge_H1)$percentage
[1] "4.47%"

compute Bayes factor

k <- bf(bridge_H1, bridge_H0)
k
Estimated Bayes factor in favor of bridge_H1 over bridge_H0: 10.42451

That is, according to Jeffreys (1961), $K>1$, in this contrast H1 is supported, then the 3-state model should be selected.

Gronau, Q. F., Sarafoglou, A., Matzke, D., Ly, A., Boehm, U., Marsman, M., ... & Steingroever, H. (2017). A tutorial on bridge sampling. Journal of mathematical psychology, 81, 80-97.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.