# Markov Chain Monte Carlo for HDP-HSMM with a Negative Binomial outcome distribution This function generates a sample from the posterior distribution of a Hidden Semi-Markov Model with a Heirarchical Dirichlet Process and a Negative Binomial outcome distribution (Johnson and Willsky, 2013). The user supplies data and priors, and a sample from the posterior distribution is returned as an mcmc object, which can be subsequently analyzed with functions provided in the coda package. ```r HDPHSMMnegbin( formula, data = parent.frame(), K = 10, b0 = 0, B0 = 1, a.alpha = 1, b.alpha = 0.1, a.gamma = 1, b.gamma = 0.1, a.omega, b.omega, e = 2, f = 2, g = 10, r = 1, burnin = 1000, mcmc = 1000, thin = 1, verbose = 0, seed = NA, beta.start = NA, P.start = NA, rho.start = NA, rho.step, nu.start = NA, omega.start = NA, gamma.start = 0.5, alpha.start = 100, ... ) ``` ## Arguments - `formula`: Model formula. - `data`: Data frame. - `K`: The number of regimes under consideration. This should be larger than the hypothesized number of regimes in the data. Note that the sampler will likely visit fewer than `K` regimes. - `b0`: The prior mean of $\beta$. This can either be a scalar or a column vector with dimension equal to the number of betas. If this takes a scalar value, then that value will serve as the prior mean for all of the betas. - `B0`: The prior precision of $\beta$. This can either be a scalar or a square matrix with dimensions equal to the number of betas. If this takes a scalar value, then that value times an identity matrix serves as the prior precision of beta. Default value of 0 is equivalent to an improper uniform prior for beta. - `a.alpha, b.alpha`: Shape and scale parameters for the Gamma distribution on $\alpha$. - `a.gamma, b.gamma`: Shape and scale parameters for the Gamma distribution on $\gamma$. - `a.omega, b.omega`: Paramaters for the Beta prior on $\omega$, which determines the regime length distribution, which is Negative Binomial, with parameters `r` and `omega`. - `e`: The hyperprior for the distribution $\rho$ See details. - `f`: The hyperprior for the distribution $\rho$. See details. - `g`: The hyperprior for the distribution $\rho$. See details. - `r`: Parameter of the Negative Binomial prior for regime durations. It is the target number of successful trials. Must be strictly positive. Higher values increase the variance of the duration distributions. - `burnin`: The number of burn-in iterations for the sampler. - `mcmc`: The number of Metropolis iterations for the sampler. - `thin`: The thinning interval used in the simulation. The number of mcmc iterations must be divisible by this value. - `verbose`: A switch which determines whether or not the progress of the sampler is printed to the screen. If `verbose` is greater than 0 the iteration number, the current beta vector, and the Metropolis acceptance rate are printed to the screen every `verbose`th iteration. - `seed`: The seed for the random number generator. If NA, the Mersenne Twister generator is used with default seed 12345; if an integer is passed it is used to seed the Mersenne twister. The user can also pass a list of length two to use the L'Ecuyer random number generator, which is suitable for parallel computation. The first element of the list is the L'Ecuyer seed, which is a vector of length six or NA (if NA a default seed of `rep(12345,6)` is used). The second element of list is a positive substream number. See the MCMCpack specification for more details. - `beta.start`: The starting value for the $\beta$ vector. This can either be a scalar or a column vector with dimension equal to the number of betas. If this takes a scalar value, then that value will serve as the starting value for all of the betas. The default value of NA will use the maximum likelihood estimate of $\beta$ as the starting value for all regimes. - `P.start`: Initial transition matrix between regimes. Should be a `K` by `K` matrix. If not provided, the default value will be uniform transition distributions. - `rho.start`: The starting value for the $\rho$ variable. This can either be a scalar or a column vector with dimension equal to the number of regimes. If the value is scalar, it will be used for all regimes. The default value is a vector of ones. - `rho.step`: Tuning parameter for the slice sampling approach to sampling $rho$. Determines the size of the step-out used to find the correct slice to draw from. Lower values are more accurate, but will take longer (up to a fixed searching limit). Default is 0.1. - `nu.start`: The starting values for the random effect, $\nu$. The default value is a vector of ones. - `omega.start`: A vector of starting values for the probability of success parameter in the Negative Binomial distribution that governs the duration distributions. - `alpha.start, gamma.start`: Scalar starting values for the $\alpha$, and $\gamma$ parameters. - `...`: further arguments to be passed. ## Returns An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package. ## Details `HDPHSMMnegbin` simulates from the posterior distribution of a HDP-HSMM with a Negative Binomial outcome distribution, allowing for multiple, arbitrary changepoints in the model. The details of the model are discussed in Johnson & Willsky (2013). The implementation here is based on a weak-limit approximation, where there is a large, though finite number of regimes that can be switched between. Unlike other changepoint models in `MCMCpack`, the HDP-HSMM approach allows for the state sequence to return to previous visited states. The model takes the following form, where we show the fixed-limit version: $$ y_t \sim \mathcal{P}oisson(\nu_t\mu_t) $$ $$ \mu_t = x_t ' \beta_k,\;\; k = 1, \ldots, K $$ $$ \nu_t \sim \mathcal{G}amma(\rho_k, \rho_k) $$ Where $K$ is an upper bound on the number of states and $\beta_k$ and $\rho_k$ are parameters when a state is $k$ at $t$. In the HDP-HSMM, there is a super-state sequence that, for a given observation, is drawn from the transition distribution and then a duration is drawn from a duration distribution to determin how long that state will stay active. After that duration, a new super-state is drawn from the transition distribution, where self-transitions are disallowed. The transition probabilities between states are assumed to follow a heirarchical Dirichlet process: $$ \pi_k \sim \mathcal{D}irichlet(\alpha\delta_1, \ldots ,\alpha\delta_K) $$ $$ \delta \sim \mathcal{D}irichlet(\gamma/K, \ldots, \gamma/K) $$ In the algorithm itself, these $\pi$ vectors are modified to remove self-transitions as discussed above. There is a unique duration distribution for each regime with the following parameters: $$ D_k \sim \mathcal{N}egBin(r, \omega_k) $$ $$ \omega_k \sim \mathcal{B}eta(a_{\omega,k}, b_{\omega, k}) $$ We assume Gaussian distribution for prior of $\beta$: $$ \beta_k \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, K $$ The overdispersion parameters have a prior with the following form: $$ f(\rho_k|e,f,g) \propto \rho^{e-1}(\rho + g)^{-(e+f)} $$ The model is simulated via blocked Gibbs conditonal on the states. The $\beta$ being simulated via the auxiliary mixture sampling method of Fuerhwirth-Schanetter et al. (2009). The $\rho$ is updated via slice sampling. The $\nu_t$ are updated their (conjugate) full conditional, which is also Gamma. The states and their durations are drawn as in Johnson & Willsky (2013). ## Examples ```r ## Not run: n <- 150 reg <- 3 true.s <- gl(reg, n/reg, n) rho.true <- c(1.5, 0.5, 3) b1.true <- c(1, -2, 2) x1 <- runif(n, 0, 2) nu.true <- rgamma(n, rho.true[true.s], rho.true[true.s]) mu <- nu.true * exp(1 + x1 * b1.true[true.s]) y <- rpois(n, mu) posterior <- HDPHSMMnegbin(y ~ x1, K = 10, verbose = 1000, e = 2, f = 2, g = 10, b0 = 0, B0 = 1/9, a.omega = 1, b.omega = 100, r = 1, rho.step = rep(0.75, times = 10), seed = list(NA, 2), omega.start = 0.05, gamma.start = 10, alpha.start = 5) plotHDPChangepoint(posterior, ylab="Density", start=1) ## End(Not run) ``` ## References Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', **Journal of Statistical Software**. 42(9): 1-21. tools:::Rd_expr_doi("10.18637/jss.v042.i09") . Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. **Scythe Statistical Library 1.0.** [http://scythe.wustl.edu.s3-website-us-east-1.amazonaws.com/](http://scythe.wustl.edu.s3-website-us-east-1.amazonaws.com/). Sylvia Fruehwirth-Schnatter, Rudolf Fruehwirth, Leonhard Held, and Havard Rue. 2009. ``Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data'', **Statistics and Computing** 19(4): 479-492. Matthew Blackwell. 2017. ``Game Changers: Detecting Shifts in Overdispersed Count Data,'' **Political Analysis** 26(2), 230-239. Matthew J. Johnson and Alan S. Willsky. 2013. ``Bayesian Nonparametric Hidden Semi-Markov Models.'' **Journal of Machine Learning Research**, 14(Feb), 673-701. ## See Also `MCMCnegbinChange`, `HDPHMMnegbin`,

HDPHSMMnegbin function