# Bayes Factor Poisson-Hidden Markov Model

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)

stan.data <- list(n=dim(dat), 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 15.068 0.053 0.953 13.121 14.497 15.184 15.736 16.641 324 0.998 lambda 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 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. • If there are some excerpts from the book, they should be clearly indicated as quotes. – Xi'an Feb 11 '19 at 9:06 • 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. – Xi'an Feb 11 '19 at 9:42 ## 2 Answers 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. 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), m=2, x=dat$V2),
iter = 1000, chains = 1)

model_H1 <- stan(model_code = PHMM,
data = list(n=dim(dat), 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.44%" error_measures(bridge_H1)$$percentage
 "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.