<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>An Introduction to BayesHMM • BayesHMM</title>
<!-- jquery --><script src="https://code.jquery.com/jquery-3.1.0.min.js" integrity="sha384-nrOSfDHtoPMzJHjVTdCopGqIqeYETSXhZDFyniQ8ZHcVy08QesyHcnOUpMpqnmWq" crossorigin="anonymous"></script><!-- Bootstrap --><link href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-BVYiiSIFeK1dGmJRAkycuHAHRg32OmUcww7on3RYdg4Va+PmSTsz/K68vbdEjh4u" crossorigin="anonymous">
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js" integrity="sha384-Tc5IQib027qvyjSMfHjOMaLkfuWVxZxUPnCJA7l2mCWNIpG9mGCD8wGNIcPD7Txa" crossorigin="anonymous"></script><!-- Font Awesome icons --><link href="https://maxcdn.bootstrapcdn.com/font-awesome/4.6.3/css/font-awesome.min.css" rel="stylesheet" integrity="sha384-T8Gy5hrqNKT+hzMclPo118YTQO6cYprQmhrYwIiQ/3axmI1hQomh7Ud2hPOy8SP1" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script><!-- sticky kit --><script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><meta property="og:title" content="An Introduction to BayesHMM">
<meta property="og:description" content="">
<meta name="twitter:card" content="summary">
<!-- mathjax --><script src="https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body>
    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">BayesHMM</a>
        <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.0.1</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../index.html">
    <span class="fa fa-home fa-lg"></span>
     
  </a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
    Articles
     
    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/introduction.html">An Introduction to BayesHMM</a>
    </li>
    <li>
      <a href="../articles/introduction2.html">An Introduction to BayesHMM</a>
    </li>
  </ul>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right">
<li>
  <a href="https://github.com/luisdamiano/BayesHMM/">
    <span class="fa fa-github fa-lg"></span>
     
  </a>
</li>
      </ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->

      
      </header><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1>An Introduction to BayesHMM</h1>
                        <h4 class="author">Luis Damiano, Michael Weylandt, Brian Peterson</h4>
            
            <h4 class="date">2019-01-03</h4>
      
      <small class="dont-index">Source: <a href="https://github.com/luisdamiano/BayesHMM//blob/master/vignettes/introduction.Rmd"><code>vignettes/introduction.Rmd</code></a></small>
      <div class="hidden name"><code>introduction.Rmd</code></div>

    </div>

    
    
<p>BayesHMM is an R Package to run full Bayesian inference on Hidden Markov Models (HMM) using the probabilistic programming language Stan. By providing an intuitive, expressive yet flexible input interface, we enable researchers to profit the most out of the Bayesian workflow. We provide the user with an expressive interface to mix and match a wide array of options for the observation and latent models, including ample choices of densities, priors, and link functions whenever covariates are present. The software enables users to fit HMM with time-homogeneous transitions as well as time-varying transition probabilities. Priors can be set for every model parameter. Implemented inference algorithms include forward (filtering), forward-backwards (smoothing), Viterbi (most likely hidden path), prior predictive sampling, and posterior predictive sampling. Graphs, tables and other convenience methods for convergence diagnosis, goodness of fit, and data analysis are provided.</p>
<p>This vignette introduces the software and briefly reviews current capabilities. Although we walk the reader through most of the functionalities, we do not discuss function arguments, implementation details, and other details typically reviewed in the documentations. A list of future features may be found in the <a href="https://github.com/luisdamiano/BayesHMM/wiki/Roadmap">Roadmap wiki page</a>.</p>
<p><strong>NOTE</strong>: Our software is work in progress. We will review naming conventions as well as other design decisions at a later stage. Fully detailed documentation will be available at the time of release.</p>
<div id="introduction" class="section level1">
<h1 class="hasAnchor">
<a href="#introduction" class="anchor"></a>Introduction</h1>
</div>
<div id="hidden-markov-models" class="section level1">
<h1 class="hasAnchor">
<a href="#hidden-markov-models" class="anchor"></a>Hidden Markov Models</h1>
<p>Real-world processes produce observable outputs characterized as signals. These can be discrete or continuous in nature, be pure or contaminated with noise, come from a stationary or non stationary source, among many other variations. These signals are modelled to allow for both theoretical descriptions and practical applications. The model itself can be deterministic or stochastic, in which case the signal is well characterized as a parametric random process whose parameters can be estimated in a well-defined manner.</p>
<p>Autocorrelation, a key feature in most signals, can be modelled in countless forms. While certainly pertinent to this purpose, high-order Markov chains can prove inconvenient when the range of the correlation amongst the observations is long. A more parsimonious approach assumes that the observed sequence is a noisy observation of an underlying hidden process represented as a first-order Markov chain. In other terms, long-range dependencies between observations are mediated via latent variables. It is important to note that the Markov property is only assumed for the hidden states and not for the observations themselves.</p>
<p>Hidden Markov Models (HMM) involve two interconnected models. The state model consists of a discrete-time, discrete-state first-order Markov chain <span class="math inline">\(z_t \in \{1, \dots, K\}\)</span> that transitions according to <span class="math inline">\(p(z_t | z_{t-1})\)</span>. In turns, the observation model is governed by <span class="math inline">\(p({\mathbf{y}}_t | z_t)\)</span>, where <span class="math inline">\({\mathbf{y}}_t \in {\mathbb{R}}^R\)</span> are the observations, emissions or output. The corresponding joint distribution is</p>
<p><span class="math display">\[
p({\mathbf{z}}_{1:T}, {\mathbf{y}}_{1:T})
  = p({\mathbf{z}}_{1:T}) p({\mathbf{y}}_{1:T} | {\mathbf{z}}_{1:T})
  = \left[ p(z_1) \prod_{t=2}^{T}{p(z_t | z_{t-1})} \right] \left[ \prod_{t=1}^{T}{p({\mathbf{y}}_t | z_{t})} \right].
\]</span></p>
<p>The non-stochastic quantities of the model are the length of the observed sequence <span class="math inline">\(T\)</span> and the number of hidden states <span class="math inline">\(K\)</span>. The observed sequence <span class="math inline">\({\mathbf{y}}_t\)</span> is a stochastic known quantity. The parameters of the models are <span class="math inline">\({\mathbf{\theta}} = ({\mathbf{\pi}}, {\mathbf{\theta}}_t, {\mathbf{\theta}}_o)\)</span>, where <span class="math inline">\({\mathbf{\pi}}\)</span> are the probabilities of the initial state distribution, <span class="math inline">\({\mathbf{\theta}}_t\)</span> are the parameters of the transition model and <span class="math inline">\({\mathbf{\theta}}_o\)</span> are the parameters of the state-conditional density function <span class="math inline">\(p({\mathbf{y}}_t | z_t)\)</span>. Their forms depend on the characteristics of the model specified by the user.</p>
<div id="observation-model" class="section level2">
<h2 class="hasAnchor">
<a href="#observation-model" class="anchor"></a>Observation model</h2>
<p>The observation model is specified by the density <span class="math inline">\(p({\mathbf{y}}_t | z_{t})\)</span>. When the output is discrete, it commonly takes the form of an observation matrix</p>
<p><span class="math display">\[
p({\mathbf{y}}_t | z_t = k, {\mathbf{\theta}}) = \text{Categorical}({\mathbf{y}}_t | {\mathbf{\theta}}_k).
\]</span></p>
<p>If the output is continuous, observations may follow for example a conditional Gaussian distribution</p>
<p><span class="math display">\[
p({\mathbf{y}}_t | z_t = k, {\mathbf{\theta}}) = \mathcal{N}({\mathbf{y}}_t | {\mathbf{\mu}}_k, {\mathbf{\Sigma}}_k).
\]</span></p>
<p>Alternatively, time-varying covariates <span class="math inline">\({\mathbf{x}}_t \in {\mathbb{R}}^M\)</span> may be used to drive the location parameter of the chosen density,</p>
<p><span class="math display">\[
p({\mathbf{y}}_t | {\mathbf{x}}_t, z_t = k, {\mathbf{\theta}}) = \mathcal{N}({\mathbf{y}}_t | {\mathbf{x}}_t {\mathbf{\beta}}^x_k, {\mathbf{\Sigma}}_k).
\]</span></p>
</div>
<div id="transition-model" class="section level2">
<h2 class="hasAnchor">
<a href="#transition-model" class="anchor"></a>Transition model</h2>
<p>In the most common case of time-homogeneous HMMs, state transitions are characterized by the <span class="math inline">\(K \times K\)</span> sized transition matrix with simplex rows <span class="math inline">\({\mathbf{A}} = \{a_{ij}\}\)</span> with <span class="math inline">\(a_{ij} = p(z_t = j | z_{t-1} = i)\)</span>. There are <span class="math inline">\(K \times (K-1)\)</span> free parameters as the rows of the matrix must sum to one.</p>
<p>Alternatively, time-varying covariates <span class="math inline">\({\mathbf{u}}_t \in {\mathbb{R}}^P\)</span> may be used to drive the location parameter of the chosen density, effectively allowing for time-varying transition probabilities. The model involves a multinomial regression whose parameters depend on the previous state taking the value <span class="math inline">\(i\)</span>. Using the softmax transform, for example,</p>
<p><span class="math display">\[
p(z_t | {\mathbf{u}}_{t}, z_{t-1} = i) = \text{softmax}({\mathbf{u}}_{t} {\mathbf{\beta}}^u_i).
\]</span></p>
</div>
<div id="initial-distribution-model" class="section level2">
<h2 class="hasAnchor">
<a href="#initial-distribution-model" class="anchor"></a>Initial distribution model</h2>
<p>Most frequently in applications, the initial distribution model is characterized by a <span class="math inline">\(K\)</span> sized simplex <span class="math inline">\({\mathbf{\pi}}\)</span> with initial probabilities <span class="math inline">\(\pi_{i} = p(z_1 = i)\)</span>. Alternatively, the first observation may be assigned to a state via a multinomial regression if relevant information is available in the form of a covariate vector <span class="math inline">\({\mathbf{v}} \in {\mathbb{R}}^Q\)</span>,</p>
<p><span class="math display">\[
p(z_1 | {\mathbf{v}}, {\mathbf{\theta}}) = \text{softmax}({\mathbf{v}} {\mathbf{\beta}}^v).
\]</span></p>
</div>
<div id="inference" class="section level2">
<h2 class="hasAnchor">
<a href="#inference" class="anchor"></a>Inference</h2>
<p>There are several quantities of interest that can be inferred via different algorithms. Our code contains the implementation of the most relevant methods for unsupervised data: forward, forward-backward and Viterbi decoding algorithms. We acknowledge the authors of the Stan Manual for the thorough illustrations and code snippets, some of which served as a starting point for our own code. As estimation is treated later, we assume that model parameters <span class="math inline">\({\mathbf{\theta}}\)</span> are known.</p>
<table class="table">
<caption>Summary of the hidden quantities and their corresponding inference algorithm. <span class="math inline">\(\dagger\)</span> Time complexity is reduced to <span class="math inline">\(O(KT)\)</span> for inference on a left-to-right (upper triangular) transition matrix.</caption>
<thead><tr class="header">
<th align="left">Name</th>
<th align="left">Hidden Quantity</th>
<th align="left">Availability at</th>
<th align="left">Algorithm</th>
<th align="left">Complexity</th>
</tr></thead>
<tbody>
<tr class="odd">
<td align="left">Filtering</td>
<td align="left"><span class="math inline">\(p(z_t | {\mathbf{y}}_{1:t})\)</span></td>
<td align="left">
<span class="math inline">\(t\)</span> (online)</td>
<td align="left">Forward</td>
<td align="left">
<span class="math inline">\(O(K^2T)\)</span>  <span class="math inline">\(O(KT)\dagger\)</span>
</td>
</tr>
<tr class="even">
<td align="left">Smoothing</td>
<td align="left"><span class="math inline">\(p(z_t | {\mathbf{y}}_{1:T})\)</span></td>
<td align="left">
<span class="math inline">\(T\)</span> (offline)</td>
<td align="left">Forward-backward</td>
<td align="left">
<span class="math inline">\(O(K^2T)\)</span>  <span class="math inline">\(O(KT)\dagger\)</span>
</td>
</tr>
<tr class="odd">
<td align="left">Fixed lag smoothing</td>
<td align="left">
<span class="math inline">\(p(z_{t - \ell} | {\mathbf{y}}_{1:t})\)</span>, <span class="math inline">\(\ell \ge 1\)</span>
</td>
<td align="left">
<span class="math inline">\(t + \ell\)</span> (lagged)</td>
<td align="left">Forward-backward</td>
<td align="left">
<span class="math inline">\(O(K^2T)\)</span>  <span class="math inline">\(O(KT)\dagger\)</span>
</td>
</tr>
<tr class="even">
<td align="left">State prediction</td>
<td align="left">
<span class="math inline">\(p(z_{t+h} | {\mathbf{y}}_{1:t})\)</span>, <span class="math inline">\(h\ge 1\)</span>
</td>
<td align="left"><span class="math inline">\(t\)</span></td>
<td align="left"></td>
<td align="left"></td>
</tr>
<tr class="odd">
<td align="left">Observation prediction</td>
<td align="left">
<span class="math inline">\(p(y_{t+h} | {\mathbf{y}}_{1:t})\)</span>, <span class="math inline">\(h\ge 1\)</span>
</td>
<td align="left"><span class="math inline">\(t\)</span></td>
<td align="left"></td>
<td align="left"></td>
</tr>
<tr class="even">
<td align="left">MAP Estimation</td>
<td align="left"><span class="math inline">\({\arg\!\max}_{{\mathbf{z}}_{1:T}} p({\mathbf{z}}_{1:T} | {\mathbf{y}}_{1:T})\)</span></td>
<td align="left"><span class="math inline">\(T\)</span></td>
<td align="left">Viterbi decoding</td>
<td align="left"><span class="math inline">\(O(K^2T)\)</span></td>
</tr>
<tr class="odd">
<td align="left">Log likelihood</td>
<td align="left"><span class="math inline">\(p({\mathbf{y}}_{1:T})\)</span></td>
<td align="left"><span class="math inline">\(T\)</span></td>
<td align="left">Forward</td>
<td align="left">
<span class="math inline">\(O(K^2T)\)</span>  <span class="math inline">\(O(KT)\dagger\)</span>
</td>
</tr>
</tbody>
</table>
<div id="filtering" class="section level3">
<h3 class="hasAnchor">
<a href="#filtering" class="anchor"></a>Filtering</h3>
<p>A filter infers the posterior distribution of the hidden states at a given step <span class="math inline">\(t\)</span> based on all the information available up to that point <span class="math inline">\(p(z_t | {\mathbf{y}}_{1:t})\)</span>. It achieves better noise reduction than simply estimating the hidden state based on the current estimate <span class="math inline">\(p(z_t | {\mathbf{y}}_{t})\)</span>. The filtering process can be run online, or recursively, as new data streams in.</p>
<p>Filtered marginals can be computed recursively using the forward algorithm <span class="citation">[@baum1967inequality]</span>. Let <span class="math inline">\(\psi_t(j) = p({\mathbf{y}}_t | z_t = j)\)</span> be the local evidence at step <span class="math inline">\(t\)</span> and <span class="math inline">\(\Psi(i, j) = p(z_t = j | z_{t-1} = i)\)</span> be the transition probability. First, the one-step-ahead predictive density is computed</p>
<p><span class="math display">\[
p(z_t = j | {\mathbf{y}}_{1:t-1}) = \sum_{i}{\Psi(i, j) p(z_{t-1} = i | {\mathbf{y}}_{1:t-1})}.
\]</span></p>
<p>Acting as prior information, this quantity is updated with observed data at the step <span class="math inline">\(t\)</span> using the Bayes rule,</p>
<span class="math display">\[\begin{align*}
\label{eq:filtered-beliefstate}
\alpha_t(j) 
  &amp; \triangleq  p(z_t = j | \mat{y}_{1:t}) \\
  &amp;= p(z_t = j | \mat{y}_{t}, \mat{y}_{1:t-1}) \\
  &amp;= Z_t^{-1} \psi_t(j) p(z_t = j | \mat{y}_{1:t-1}) \\
\end{align*}\]</span>
<p>where the normalization constant is given by</p>
<p><span class="math display">\[
Z_t
  \triangleq  p({\mathbf{y}}_t | {\mathbf{y}}_{1:t-1})
  = \sum_{l=1}^{K}{p({\mathbf{y}}_{t} | z_t = l) p(z_t = l | {\mathbf{y}}_{1:t-1})}
  = \sum_{l=1}^{K}{\psi_t(l) p(z_t = l | {\mathbf{y}}_{1:t-1})}.
\]</span></p>
<p>This predict-update cycle results in the filtered belief states at step <span class="math inline">\(t\)</span>. As this algorithm only requires the evaluation of the quantities <span class="math inline">\(\psi_t(j)\)</span> for each value of <span class="math inline">\(z_t\)</span> for every <span class="math inline">\(t\)</span> and fixed <span class="math inline">\({\mathbf{y}}_t\)</span>, the posterior distribution of the latent states is independent of the form of the observation density or indeed of whether the observed variables are continuous or discrete <span class="citation">[@jordan2003introduction]</span>. In other words, <span class="math inline">\(\alpha_t(j)\)</span> does not depend on the complete form of the density function <span class="math inline">\(p({\mathbf{y}} | {\mathbf{z}})\)</span> but only on the point values <span class="math inline">\(p({\mathbf{y}}_t | z_t = j)\)</span> for every <span class="math inline">\({\mathbf{y}}_t\)</span> and <span class="math inline">\(z_t\)</span>.</p>
<p>Let <span class="math inline">\({\mathbf{\alpha}}_t\)</span> be a <span class="math inline">\(K\)</span>-sized vector with the filtered belief states at step <span class="math inline">\(t\)</span>, <span class="math inline">\({\mathbf{\psi}}_t(j)\)</span> be the <span class="math inline">\(K\)</span>-sized vector of local evidence at step <span class="math inline">\(t\)</span>, <span class="math inline">\({\mathbf{\Psi}}\)</span> be the transition matrix and <span class="math inline">\({\mathbf{u}} \odot {\mathbf{v}}\)</span> be the Hadamard product, representing element-wise vector multiplication. Then, the Bayesian updating procedure can be expressed in matrix notation as</p>
<p><span class="math display">\[
{\mathbf{\alpha}}_t \propto {\mathbf{\psi}}_t \odot ({\mathbf{\Psi}}^T {\mathbf{\alpha}}_{t-1}).
\]</span></p>
<p>In addition to computing the hidden states, the algorithm yields the log likelihood</p>
<p><span class="math display">\[
{\mathcal{L}}= \log p({\mathbf{y}}_{1:T} | {\mathbf{\theta}}) = \sum_{t=1}^{T}{\log p({\mathbf{y}}_{t} | {\mathbf{y}}_{1:t-1})} = \sum_{t=1}^{T}{\log Z_t}.
\]</span></p>
<p>The time complexity of the algorithm is <span class="math inline">\(O(K^2T)\)</span>: there are <span class="math inline">\(K \times K\)</span> iterations within each of the <span class="math inline">\(T\)</span> iterations of the outer loop. Brute-forcing through all possible hidden states <span class="math inline">\(K^T\)</span> would prove prohibitive for realistic problems as time complexity increases exponentially with sequence length <span class="math inline">\(O(K^TT)\)</span>.</p>
</div>
<div id="smoothing" class="section level3">
<h3 class="hasAnchor">
<a href="#smoothing" class="anchor"></a>Smoothing</h3>
<p>A smoother infers the posterior distribution of the hidden states at a given state based on all the observations or evidence <span class="math inline">\(p(z_t | {\mathbf{y}}_{1:T})\)</span>. Although noise and uncertainty are significantly reduced as a result of conditioning on past and future data, the smoothing process can only be run offline.</p>
<p>Inference can be done by means of the forward-backward algorithm, which also plays an important role in the Baum-Welch algorithm for learning model parameters <span class="citation">[@baum1967inequality; @baum1970maximization]</span>. Let <span class="math inline">\(\gamma_t(j)\)</span> be the desired smoothed posterior marginal,</p>
<p><span class="math display">\[
\gamma_t(j)
  \triangleq p(z_t = j | {\mathbf{y}}_{1:T}),
\]</span></p>
<p><span class="math inline">\(\alpha_t(j)\)</span> be the filtered belief state at the step <span class="math inline">\(t\)</span> as defined previously, and <span class="math inline">\(\beta_t(j)\)</span> be the conditional likelihood of future evidence given that the hidden state at step <span class="math inline">\(t\)</span> is <span class="math inline">\(j\)</span>,</p>
<p><span class="math display">\[
\beta_t(j) 
  \triangleq p({\mathbf{y}}_{t+1:T} | z_t = j).
\]</span></p>
<p>Then, the chain of smoothed marginals can be segregated into the past and the future components by conditioning on the belief state <span class="math inline">\(z_t\)</span>,</p>
<p><span class="math display">\[
\gamma_t(j)
  = \frac{\alpha_t(j) \beta_t(j)}{p({\mathbf{y}}_{1:T})}
  \propto \alpha_t(j) \beta_t(j).
\]</span></p>
<p>The future component can be computed recursively from right to left:</p>
<span class="math display">\[\begin{align*}
\beta_{t-1}(i)
  &amp;= p(\mat{y}_{t:T} | z_{t-1} = i) \\
  &amp;= \sum_{j=1}^{K}{p(z_t =j, \mat{y}_{t}, \mat{y}_{t+1:T} | z_{t-1} = i)} \\
  &amp;= \sum_{j=1}^{K}{p(\mat{y}_{t+1:T} | z_t = j)p(z_t = j, \mat{y}_{t} | z_{t-1} = i)} \\
  &amp;= \sum_{j=1}^{K}{p(\mat{y}_{t+1:T} | z_t = j)p(\mat{y}_t | z_t = j)p(z_t = j | z_{t-1} = i)} \\
  &amp;= \sum_{j=1}^{K}{\beta_t(j) \psi_t(j) \Psi(i, j)}
\end{align*}\]</span>
<p>Let <span class="math inline">\({\mathbf{\beta}}_t\)</span> be a <span class="math inline">\(K\)</span>-sized vector with the conditional likelihood of future evidence given the hidden state at step <span class="math inline">\(t\)</span>. Then, the backward procedure can be expressed in matrix notation as</p>
<p><span class="math display">\[
{\mathbf{\beta}}_{t-1} \propto {\mathbf{\Psi}} ({\mathbf{\psi}}_t \odot {\mathbf{\beta}}_{t}).
\]</span></p>
<p>At the last step, the base case is given by <span class="math display">\[
\beta_{T}(i)
  = p({\mathbf{y}}_{T+1:T} | z_{T} = i) = p(\varnothing | z_t = i) = 1.
\]</span></p>
<p>Intuitively, the forward-backward algorithm passes information from left to right and then from right to left, combining them at each node. A straightforward implementation of the algorithm runs in <span class="math inline">\(O(K^2 T)\)</span> time because of the <span class="math inline">\(K \times K\)</span> matrix multiplication at each step. Nonetheless the frequent description as two subsequent passes, these procedures are not inherently sequential and share no information. As a result, they could be implemented in parallel.</p>
<p>There is a significant reduction if the transition matrix is sparse. Inference for a left-to-right (upper triangular) transition matrix, a model where the state index increases or stays the same as time passes, runs in <span class="math inline">\(O(TK)\)</span> time <span class="citation">[@bakis1976continuous;@jelinek1976continuous]</span>. Additional assumptions about the form of the transition matrix may ease complexity further, for example decreasing the time to <span class="math inline">\(O(TK\log K)\)</span> if <span class="math inline">\(\psi(i, j) \propto \exp(-\sigma^2 |{\mathbf{z}}_i - {\mathbf{z}}_j|)\)</span>. Finally, ad-hoc data pre-processing strategies may help control complexity, for example by pruning nodes with low conditional probability of occurrence.</p>
<p>The forward-backward algorithm was designed to exploit via recursion the conditional independencies in the HMM. First, the posterior marginal probability of the latent states at a given time step is broken down into two quantities: the past and the future components. Second, taking advantage of the Markov properties, each of the two are further broken down into simpler pieces via conditioning and marginalizing, thus creating an efficient predict-update cycle.</p>
<p>This strategy makes otherwise unfeasible calculations possible. Consider for example a time series with <span class="math inline">\(T = 100\)</span> observations and <span class="math inline">\(K = 5\)</span> hidden states. Summing the joint probability over all possible state sequences would involve <span class="math inline">\(2 \times 100 \times 5^{100} \approx 10^{72}\)</span> computations, while the forward and backward passes only take <span class="math inline">\(3,000\)</span> each. Moreover, one pass may be avoided depending on the goal of the data analysis. Summing the forward probabilities at the last time step is enough to compute the likelihood, and the backwards recursion would be only needed if the posterior probabilities of the states were also required.</p>
</div>
<div id="map-viterbi" class="section level3">
<h3 class="hasAnchor">
<a href="#map-viterbi" class="anchor"></a>MAP: Viterbi</h3>
<p>It is also of interest to compute the most probable state sequence or path,</p>
<p><span class="math display">\[
{\mathbf{z}}^* = {\arg\!\max}_{{\mathbf{z}}_{1:T}} p({\mathbf{z}}_{1:T} | {\mathbf{y}}_{1:T}).
\]</span></p>
<p>The jointly most probable sequence of states can be inferred via maximum <em>a posteriori</em> (MAP) estimation. Note that the jointly most probable sequence is not necessarily the same as the sequence of marginally most probable states given by the maximizer of the posterior marginals (MPM),</p>
<p><span class="math display">\[
{\mathbf{\hat{z}}} = ({\arg\!\max}_{z_1} p(z_1 | {\mathbf{y}}_{1:T}), \ \dots, \ {\arg\!\max}_{z_T} p(z_T | {\mathbf{y}}_{1:T})),
\]</span></p>
<p>which maximizes the expected number of correct individual states.</p>
<p>The MAP estimate is always globally consistent: while locally a state may be most probable at a given step, the Viterbi or max-sum algorithm decodes the most likely single plausible path <span class="citation">[@viterbi1967error]</span>. Furthermore, the MPM sequence may have zero joint probability if it includes two successive states that, while being individually the most probable, are connected in the transition matrix by a zero. On the other hand, MPM may be considered more robust since the state at each step is estimated by averaging over its neighbors rather than conditioning on a specific value of them.</p>
<p>The Viterbi applies the max-sum algorithm in a forward fashion plus a traceback procedure to recover the most probable path. In simple terms, once the most probable state <span class="math inline">\(z_t\)</span> is estimated, the procedure conditions the previous states on it. Let <span class="math inline">\(\delta_t(j)\)</span> be the probability of arriving to the state <span class="math inline">\(j\)</span> at step <span class="math inline">\(t\)</span> given the most probable path was taken,</p>
<p><span class="math display">\[
\delta_t(j)
  \triangleq \max_{z_1, \dots, z_{t-1}} p({\mathbf{z}}_{1:t-1}, z_t = j | {\mathbf{y}}_{1:t}).
\]</span></p>
<p>The most probable path to state <span class="math inline">\(j\)</span> at step <span class="math inline">\(t\)</span> consists of the most probable path to some other state <span class="math inline">\(i\)</span> at point <span class="math inline">\(t-1\)</span>, followed by a transition from <span class="math inline">\(i\)</span> to <span class="math inline">\(j\)</span>,</p>
<p><span class="math display">\[
\delta_t(j)
  = \max_{i} \delta_{t-1}(i) \ \psi(i, j) \ \psi_t(j).
\]</span></p>
<p>Additionally, the most likely previous state on the most probable path to <span class="math inline">\(j\)</span> at step <span class="math inline">\(t\)</span> is given by <span class="math display">\[
a_t(j)
  = {\arg\!\max}_{i} \delta_{t-1}(i) \ \psi(i, j) \ \psi_t(j).
\]</span></p>
<p>By initializing with <span class="math inline">\(\delta_1 = \pi_j \phi_1(j)\)</span> and terminating with the most probable final state <span class="math inline">\(z_T^* = {\arg\!\max}_{i} \delta_T(i)\)</span>, the most probable sequence of states is estimated using the traceback,</p>
<p><span class="math display">\[
z_t^* = a_{t+1}(z_{t+1}^*).
\]</span></p>
<p>It is advisable to work in the log domain to avoid numerical underflow,</p>
<p><span class="math display">\[
\delta_t(j)
  \triangleq \max_{{\mathbf{z}}_{1:t-1}} \log p({\mathbf{z}}_{1:t-1}, z_t = j | {\mathbf{y}}_{1:t})
  = \max_{i} \log \delta_{t-1}(i) + \log \psi(i, j) + \log \psi_t(j).
\]</span></p>
<p>As with the backward-forward algorithm, the time complexity of Viterbi is <span class="math inline">\(O(K^2T)\)</span> and the space complexity is <span class="math inline">\(O(KT)\)</span>. If the transition matrix has the form <span class="math inline">\(\psi(i, j) \propto \exp(-\sigma^2 ||{\mathbf{z}}_i - {\mathbf{z}}_j||^2)\)</span>, implementation runs in <span class="math inline">\(O(TK)\)</span> time.</p>
</div>
</div>
<div id="parameter-estimation" class="section level2">
<h2 class="hasAnchor">
<a href="#parameter-estimation" class="anchor"></a>Parameter estimation</h2>
<p>The model likelihood can be derived from the definition of the quantity <span class="math inline">\(\gamma_t(j)\)</span>: given that its sum over all possible values of the latent variable must equal one, the log likelihood at time index <span class="math inline">\(t\)</span> becomes</p>
<p><span class="math display">\[
{\mathcal{L}}_t = \sum_{i = 1}^{K}{\alpha_t(i) \beta_{t}(i)}.
\]</span></p>
<p>The last step <span class="math inline">\(T\)</span> has two convenient characteristics. First, the recurrent nature of the forward probability implies that the last iteration retains the information of all the intermediate state probabilities. Second, the base case for the backwards quantity is <span class="math inline">\(\beta_{T}(i) = 1\)</span>. Consequently, the log likelihood reduces to</p>
<p><span class="math display">\[
{\mathcal{L}}_T \propto \sum_{i = 1}^{K}{\alpha_T(i)}.
\]</span></p>
</div>
<div id="naming-convention" class="section level2">
<h2 class="hasAnchor">
<a href="#naming-convention" class="anchor"></a>Naming convention</h2>
<p>Documentation and software adopts the following naming convention. Minor modifications may be found due to syntax restrictions in R or Stan. Time suffix <code>_t</code> is optional.</p>
<table class="table">
<thead><tr class="header">
<th>Constants</th>
<th></th>
</tr></thead>
<tbody>
<tr class="odd">
<td>R</td>
<td>Observation dimension</td>
</tr>
<tr class="even">
<td>K</td>
<td>Number of hidden states</td>
</tr>
<tr class="odd">
<td>M</td>
<td>Number of covariates for the observation model</td>
</tr>
<tr class="even">
<td>P</td>
<td>Number of covariates for the transition model</td>
</tr>
<tr class="odd">
<td>Q</td>
<td>Number of covariates for the initial model</td>
</tr>
<tr class="even">
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>Covariates</td>
<td></td>
</tr>
<tr class="even">
<td>x_t</td>
<td>Time-varying covariates for the observation model</td>
</tr>
<tr class="odd">
<td>u_t</td>
<td>Time-varying covariates for the transition model</td>
</tr>
<tr class="even">
<td>v</td>
<td>Covariates for the initial model</td>
</tr>
<tr class="odd">
<td></td>
<td></td>
</tr>
<tr class="even">
<td>Known-stochastic quantities</td>
<td></td>
</tr>
<tr class="odd">
<td>y_t</td>
<td>Observation vector</td>
</tr>
<tr class="even">
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>Model parameters</td>
<td></td>
</tr>
<tr class="even">
<td>*_kr</td>
<td>Ex. mu_11 or sigma_11 (suffixes k and r are optional)</td>
</tr>
<tr class="odd">
<td>A</td>
<td>Transition model parameters (if no covariates)</td>
</tr>
<tr class="even">
<td>pi</td>
<td>Initial distribution parameters (if no covariates)</td>
</tr>
<tr class="odd">
<td>xBeta</td>
<td>Regression parameters for the observation model</td>
</tr>
<tr class="even">
<td>uBeta</td>
<td>Regression parameters for the transition model</td>
</tr>
<tr class="odd">
<td>vBeta</td>
<td>Regression parameters for the initial model</td>
</tr>
<tr class="even">
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>Estimated quantities</td>
<td></td>
</tr>
<tr class="even">
<td>z_t</td>
<td>Hidden state</td>
</tr>
<tr class="odd">
<td>alpha_t Filtered probability</td>
<td></td>
</tr>
<tr class="even">
<td>gamma_t Smoothed probability</td>
<td></td>
</tr>
<tr class="odd">
<td>zstar</td>
<td>Viterbi</td>
</tr>
<tr class="even">
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>(Prior/Posterior) predictive quantities</td>
<td></td>
</tr>
<tr class="even">
<td>yPred</td>
<td>Sample of observations drawn from the predictive density</td>
</tr>
<tr class="odd">
<td>zPred</td>
<td>Sample of latent path drawn from the predictive density</td>
</tr>
</tbody>
</table>
</div>
</div>
<div id="using-bayeshmm" class="section level1">
<h1 class="hasAnchor">
<a href="#using-bayeshmm" class="anchor"></a>Using BayesHMM</h1>
<p>The typical data analysis workflow includes the following steps: (1) specify, (2) validate, (3) simulate, (4) fit, (5) diagnose, (6) visualize, (7) compare. All steps except for 5 and 7 are implemented as of today.</p>
<p>You may install the lastest version of our software from GitHub. Please note that BayesHMM requires rstan (&gt;= 2.17), the R interface for the probabilistic programming language <a href="http://mc-stan.org/">Stan</a>.</p>
<pre><code>#&gt; Loading required package: Rcpp
#&gt; Hola! BayesHHM v0.0.1 here o/</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">devtools<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/devtools/topics/reexports">install_github</a></span>(<span class="st">"luisdamiano/BayesHMM"</span>, <span class="dt">ref =</span> <span class="st">"master"</span>)
<span class="kw">library</span>(BayesHMM)</code></pre></div>
<div id="specify" class="section level2">
<h2 class="hasAnchor">
<a href="#specify" class="anchor"></a>1. Specify</h2>
<p>A HMM may be specified by calling the <code>hmm</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">2</span>,
  <span class="dt">observation =</span> {...}
  <span class="dt">initial     =</span> {...}
  <span class="dt">transition  =</span> {...}
  <span class="dt">name =</span> <span class="st">"Model name"</span>
)</code></pre></div>
<p>where <span class="math inline">\(K\)</span> is the number of hidden states and <span class="math inline">\(R\)</span> is the number of dimensions in the observation variable.</p>
<p>The <code>observation</code>, <code>initial</code>, and <code>transition</code> fields rely on S3 objects called <code>Density</code> to specify density form, parameter priors, and fixed values for parameters. User may specify bounds in the parameter space as well as truncation in prior densities. For example:</p>
<ul>
<li>
<code><a href="../reference/Gaussian.html">Gaussian(mu = 0, sigma = 1)</a></code> specifies a Gaussian density with fixed parameters.</li>
<li>
<code><a href="../reference/Gaussian.html">Gaussian(mu = Gaussian(0, 10), sigma = Cauchy(0, 10, bounds = list(0, NA)))</a></code> specifies a Gaussian density with a Gaussian prior on the location parameter and a Cauchy prior on the scale parameter, which is bounded on <span class="math inline">\([0, \infty)\)</span>.</li>
</ul>
<p>Currently available densities are listed below:</p>
<ul>
<li>Observation model:
<ul>
<li>Univariate: Bernoulli, Beta, Binomial, Categorical, Cauchy, Dirichlet, Gaussian, Multinomial, Negative Binomial (traditional and location parameterizations), Poisson, Student.</li>
<li>Multivariate: Multivariate Gaussian (traditional, Cholesky decomposition of covariance matrix, and Cholesky decomposition of correlation matrix parameterizations), Multivariate Student.</li>
<li>Regressions: Bernoulli regression with logit link, Binomial regression (logit and probit links), Softmax regression, Gaussian regression.</li>
<li>Prior-only density: LKJ, Wishart.</li>
</ul>
</li>
<li>Transition model: Dirichlet, Softmax regression.</li>
<li>Initial model: Dirichlet, Softmax regression.</li>
</ul>
<div id="observation-model-1" class="section level3">
<h3 class="hasAnchor">
<a href="#observation-model-1" class="anchor"></a>Observation model</h3>
<p>Internally, a specification stores either <span class="math inline">\(K\)</span> multivariate densities (i.e. one multivariate density for state) or <span class="math inline">\(K \times R\)</span> univariate densities (i.e. one univariate density for each dimension in the observation variable and each state). As shown in the next table, densities are recycled to simplify user-input code.</p>
<table class="table">
<thead><tr class="header">
<th>R</th>
<th>User input</th>
<th>Density</th>
<th>Action</th>
</tr></thead>
<tbody>
<tr class="odd">
<td>1</td>
<td>1</td>
<td>UnivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="even">
<td>1</td>
<td>K</td>
<td>UnivariateDensity</td>
<td>repeat_none</td>
</tr>
<tr class="odd">
<td>&gt;1</td>
<td>1</td>
<td>MultivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="even">
<td>&gt;1</td>
<td>1</td>
<td>UnivariateDensity</td>
<td>repeat_KxR</td>
</tr>
<tr class="odd">
<td>&gt;1</td>
<td>K</td>
<td>MultivariateDensity</td>
<td>repeat_none</td>
</tr>
<tr class="even">
<td>&gt;1</td>
<td>K</td>
<td>UnivariateDensity</td>
<td>repeat_R</td>
</tr>
<tr class="odd">
<td>&gt;1</td>
<td>R</td>
<td>UnivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="even">
<td>&gt;1</td>
<td>KxR</td>
<td>UnivariateDensity</td>
<td>repeat_none</td>
</tr>
</tbody>
</table>
<p>For a univariate model (<span class="math inline">\(R = 1\)</span>):</p>
<ul>
<li>User enters one univariate density: specify the same density in each hidden state variable.</li>
<li>User enters <span class="math inline">\(K\)</span> univariate density: specify one different density in each hidden state (ex. Gaussian in one state and Student in the other).</li>
</ul>
<p>This system is flexible enough to let the user specify different densities and priors for each state. The following example demonstrates a (rather unrealistic) model where states are a priori believed to have negative, near zero, and positive location parameters.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">1</span>,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="op">-</span><span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">10</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>   <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">10</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>  <span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">10</span>),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>()),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>()),
  <span class="dt">name =</span> <span class="st">"Univariate Gaussian"</span>
)</code></pre></div>
<p>For a multivariate model (<span class="math inline">\(R &gt; 1\)</span>):</p>
<ul>
<li>User enters one univariate density: same density for each dimension of the observation variable in all states (i.e. same density and priors for the R variables in the K latent states).</li>
<li>User enters one multivariate density: same density for the observation vector in all hidden states.</li>
<li>User enters <span class="math inline">\(K\)</span> univariate density: each dimension of the observation variable has the same specification within a each state.</li>
<li>User enters <span class="math inline">\(K\)</span> multivariate density: specify a multivariate density for the observation vector for each state.</li>
<li>User enters <span class="math inline">\(R\)</span> univariate density: specify one density for each element of the observation vector, which is the same across state.</li>
<li>User enters <span class="math inline">\(R \times K\)</span> univariate density: specify density for each element of the observation vector in each state. Order is <span class="math inline">\(k_1 r_1, \dots, k_1, r_R, k_2r_1, \dots, k_2r_R, \dots, k_Kr_1, \dots, k_K,r_R\)</span>.</li>
<li>When <span class="math inline">\(K = R\)</span>, <span class="math inline">\(K\)</span> prevails.</li>
</ul>
</div>
<div id="transition-model-1" class="section level3">
<h3 class="hasAnchor">
<a href="#transition-model-1" class="anchor"></a>Transition model</h3>
<p>The user may enter a <code>Density</code> S3 object such as a <code>Dirichlet</code>, or a <code>LinkDensity</code> S3 object (inherits from <code>Density</code>) such as <code><a href="../reference/TransitionSoftmax.html">TransitionSoftmax()</a></code>. The latter allows to specify priors for the transition regression parameter vector. For example, the following example sets regularizing priors near zero on the regression coefficients (because softmax coefficients are location invariant, we use priors to better inform the model). Note that this model is hard to identify, although we achieved satisfactory results optimizing the posterior density (i.e. maximum a posterior estimates).</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">2</span>, <span class="dt">R =</span> <span class="dv">1</span>,
  <span class="dt">observation =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
    <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
    <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Student.html">Student</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">10</span>, <span class="dt">nu =</span> <span class="dv">1</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
  ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/TransitionSoftmax.html">TransitionSoftmax</a></span>(
    <span class="dt">uBeta =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">5</span>)
  ),
  <span class="dt">name =</span> <span class="st">"TVHMM Softmax Univariate Gaussian"</span>
)</code></pre></div>
<p>Internally, a specification stores either <span class="math inline">\(K\)</span> multivariate densities (i.e. one multivariate density for each row in the transition matrix) or <span class="math inline">\(K \times K\)</span> univariate densities (i.e. one univariate density for each element in the transition matrix). Densities are recycled accordingly.</p>
<table class="table">
<thead><tr class="header">
<th>User-input</th>
<th>Density</th>
<th>Action</th>
</tr></thead>
<tbody>
<tr class="odd">
<td>1</td>
<td>UnivariateDensity</td>
<td>repeat_KxK</td>
</tr>
<tr class="even">
<td>1</td>
<td>MultivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="odd">
<td>1</td>
<td>LinkDensity</td>
<td>nothing</td>
</tr>
<tr class="even">
<td>K</td>
<td>UnivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="odd">
<td>K</td>
<td>MultivariateDensity</td>
<td>repeat_none_multivariate</td>
</tr>
<tr class="even">
<td>KxK</td>
<td>UnivariateDensity</td>
<td>repeat_none_univariate</td>
</tr>
</tbody>
</table>
<p><strong>NOTE</strong>: As of the time of this writing, fixed values for transition parameters have not been implemented yet.</p>
</div>
<div id="initial-model" class="section level3">
<h3 class="hasAnchor">
<a href="#initial-model" class="anchor"></a>Initial model</h3>
<p>Internally, a specification stores either one multivariate densities for the whole initial probability vector, or <span class="math inline">\(K\)</span> univariate densities for each element of said vector. Densities are recycled accordingly.</p>
<table class="table">
<thead><tr class="header">
<th>User-input</th>
<th>Density</th>
<th>Action</th>
</tr></thead>
<tbody>
<tr class="odd">
<td>1</td>
<td>UnivariateDensity</td>
<td>repeat_K</td>
</tr>
<tr class="even">
<td>1</td>
<td>MultivariateDensity</td>
<td>repeat_none_multivariate</td>
</tr>
<tr class="odd">
<td>1</td>
<td>LinkDensity</td>
<td>repeat_none_multivariate</td>
</tr>
<tr class="even">
<td>K</td>
<td>UnivariateDensity</td>
<td>repeat_none_univariate</td>
</tr>
</tbody>
</table>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">2</span>, <span class="dt">R =</span> <span class="dv">1</span>,
  <span class="dt">observation =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
    <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
    <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Student.html">Student</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">10</span>, <span class="dt">nu =</span> <span class="dv">1</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
  ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/InitialSoftmax.html">InitialSoftmax</a></span>(
    <span class="dt">vBeta =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">5</span>)
  ),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>()),
  <span class="dt">name =</span> <span class="st">"Univariate Gaussian with covariates for initial model"</span>
)</code></pre></div>
</div>
<div id="examples" class="section level3">
<h3 class="hasAnchor">
<a href="#examples" class="anchor"></a>Examples</h3>
<p>Many possible configurations are shown in the Appendix .</p>
</div>
</div>
<div id="validate" class="section level2">
<h2 class="hasAnchor">
<a href="#validate" class="anchor"></a>2. Validate</h2>
<p>Bayesian Software Validation allows users to test software accuracy [CITE]. We provide a one-liner to run a validation protocol based on the prior predictive density. Iterations are run in parallel using <code>foreach</code> and <code>doParallel</code>.</p>
<p><strong>Validation protocol inspired by Simulation Based Calibration (cite)</strong></p>
<ul>
<li>Compile the prior predictive model (i.e. no likelihood statement in the Stan code).</li>
<li>Draw <span class="math inline">\(N\)</span> samples of the parameter vector <span class="math inline">\(\theta\)</span> and the observation vector <span class="math inline">\({\mathbf{y}}_t\)</span> from prior predictive density.</li>
<li>Compile the posterior predictive model (i.e. Stan code includes both prior density and likelihood statement).</li>
<li>For all <span class="math inline">\(n \in 1, \dots, N\)</span>:
<ul>
<li>Feed <span class="math inline">\({\mathbf{y}}_t^{(n)}\)</span> to the full model.</li>
<li>Draw one posterior sample of the observation variable <span class="math inline">\({\mathbf{y_t}}^{(n)}_{\text{new}}\)</span>.</li>
<li>Collect Hamiltonian Monte Carlo diagnostics: number of divergences, number of times max tree depth is reached, maximum leapfrogs, warm up and sample times.</li>
<li>Collect posterior sampling diagnostics: posterior summary measures (mean, sd, quantiles), comparison against the true value (rank), MCMC convergence measures (Monte Carlo SE, ESS, R Hat).</li>
<li>Collect posterior predictive diagnostics: observation ranks, Kolmogorov-Smirnov statistic for observed sample vs posterior predictive samples.</li>
<li>Other: we expect to enrich the protocol with more relevant quantities in future iterations.</li>
</ul>
</li>
</ul>
<p>There is no such a thing as an universal, automatic way to validate models. The software computes and summarizes different diagnostic measures and provides the user with high-level tools to assess calibration (including printouts, tables, and plots). Because the user may want to single out and inspect one of the iterations, in future versions it will be able to recover (or reproduce) the full posterior sample drawn in the <span class="math inline">\(n\)</span>-th iteration.</p>
<p>In the following example, we specify a HMM and validate its calibration.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">  mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
    <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">2</span>,
    <span class="dt">observation =</span>
      <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="op">-</span><span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">      </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>   <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">      </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>  <span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>),
    <span class="dt">initial     =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
    <span class="dt">transition  =</span>
      <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">1.0</span>, <span class="fl">0.2</span>, <span class="fl">0.2</span>)) <span class="op">+</span>
<span class="st">      </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.2</span>, <span class="fl">1.0</span>, <span class="fl">0.2</span>)) <span class="op">+</span>
<span class="st">      </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.2</span>, <span class="fl">0.2</span>, <span class="fl">1.0</span>)),
    <span class="dt">name =</span> <span class="st">"Dummy Model"</span>
  )

  val &lt;-<span class="st"> </span><span class="kw"><a href="../reference/validate_calibration.html">validate_calibration</a></span>(mySpec, <span class="dt">N =</span> <span class="dv">2</span>, <span class="dt">T =</span> <span class="dv">100</span>, <span class="dt">iter =</span> <span class="dv">500</span>, <span class="dt">seed =</span> <span class="dv">9000</span>)
<span class="co">#&gt; hash mismatch so recompiling; make sure Stan code ends with a blank line</span>
<span class="co">#&gt; hash mismatch so recompiling; make sure Stan code ends with a blank line</span></code></pre></div>
<p>The function returns a named list with two elements, <code>chains</code> and <code>parameters</code>, which are data frames storing the aforementioned summary quantities.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">knitr<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/knitr/topics/kable">kable</a></span>(
  <span class="kw">head</span>(val<span class="op">$</span>chains), 
  <span class="dt">digits =</span> <span class="dv">2</span>, <span class="dt">caption =</span> <span class="st">"Chains diagnostics"</span>
)</code></pre></div>
<table class="table">
<caption>Chains diagnostics</caption>
<thead><tr class="header">
<th align="left">model</th>
<th align="right">chain</th>
<th align="right">divergences</th>
<th align="right">maxTreeDepth</th>
<th align="right">maxNLeapfrog</th>
<th align="right">warmup</th>
<th align="right">sample</th>
<th align="right">r1.y25Rank</th>
<th align="right">r1.yMedRank</th>
<th align="right">r1.y75Rank</th>
<th align="right">r1.yMeanRank</th>
<th align="right">r1.ksMean</th>
<th align="right">r1.ksMedian</th>
<th align="right">r2.y25Rank</th>
<th align="right">r2.yMedRank</th>
<th align="right">r2.y75Rank</th>
<th align="right">r2.yMeanRank</th>
<th align="right">r2.ksMean</th>
<th align="right">r2.ksMedian</th>
<th align="right">seed</th>
<th align="right">n</th>
</tr></thead>
<tbody>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="right">33</td>
<td align="right">0</td>
<td align="right">25</td>
<td align="right">0.76</td>
<td align="right">0.52</td>
<td align="right">0.12</td>
<td align="right">0.61</td>
<td align="right">0.82</td>
<td align="right">0.35</td>
<td align="right">0.29</td>
<td align="right">0.30</td>
<td align="right">0.90</td>
<td align="right">0.80</td>
<td align="right">0.65</td>
<td align="right">0.67</td>
<td align="right">0.22</td>
<td align="right">0.21</td>
<td align="right">9000</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">2</td>
<td align="right">0</td>
<td align="right">0</td>
<td align="right">15</td>
<td align="right">0.81</td>
<td align="right">0.55</td>
<td align="right">0.02</td>
<td align="right">0.14</td>
<td align="right">0.74</td>
<td align="right">0.17</td>
<td align="right">0.34</td>
<td align="right">0.35</td>
<td align="right">0.83</td>
<td align="right">0.80</td>
<td align="right">0.57</td>
<td align="right">0.62</td>
<td align="right">0.23</td>
<td align="right">0.22</td>
<td align="right">9000</td>
<td align="right">1</td>
</tr>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">3</td>
<td align="right">0</td>
<td align="right">0</td>
<td align="right">15</td>
<td align="right">0.84</td>
<td align="right">0.52</td>
<td align="right">0.01</td>
<td align="right">0.17</td>
<td align="right">0.41</td>
<td align="right">0.11</td>
<td align="right">0.37</td>
<td align="right">0.37</td>
<td align="right">0.80</td>
<td align="right">0.76</td>
<td align="right">0.70</td>
<td align="right">0.65</td>
<td align="right">0.22</td>
<td align="right">0.22</td>
<td align="right">9000</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">4</td>
<td align="right">0</td>
<td align="right">0</td>
<td align="right">31</td>
<td align="right">0.98</td>
<td align="right">0.54</td>
<td align="right">0.02</td>
<td align="right">0.03</td>
<td align="right">0.45</td>
<td align="right">0.00</td>
<td align="right">0.34</td>
<td align="right">0.34</td>
<td align="right">0.79</td>
<td align="right">0.78</td>
<td align="right">0.59</td>
<td align="right">0.60</td>
<td align="right">0.23</td>
<td align="right">0.22</td>
<td align="right">9000</td>
<td align="right">1</td>
</tr>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="right">0</td>
<td align="right">0</td>
<td align="right">31</td>
<td align="right">1.23</td>
<td align="right">0.54</td>
<td align="right">0.05</td>
<td align="right">0.02</td>
<td align="right">0.61</td>
<td align="right">0.03</td>
<td align="right">0.50</td>
<td align="right">0.53</td>
<td align="right">0.00</td>
<td align="right">0.13</td>
<td align="right">0.68</td>
<td align="right">0.48</td>
<td align="right">0.23</td>
<td align="right">0.22</td>
<td align="right">9000</td>
<td align="right">2</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">2</td>
<td align="right">0</td>
<td align="right">0</td>
<td align="right">15</td>
<td align="right">1.11</td>
<td align="right">0.62</td>
<td align="right">0.02</td>
<td align="right">0.00</td>
<td align="right">0.24</td>
<td align="right">0.01</td>
<td align="right">0.50</td>
<td align="right">0.52</td>
<td align="right">0.00</td>
<td align="right">0.09</td>
<td align="right">0.82</td>
<td align="right">0.53</td>
<td align="right">0.23</td>
<td align="right">0.22</td>
<td align="right">9000</td>
<td align="right">2</td>
</tr>
</tbody>
</table>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
knitr<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/knitr/topics/kable">kable</a></span>(
  <span class="kw">head</span>(val<span class="op">$</span>parameters), 
  <span class="dt">digits =</span> <span class="dv">2</span>, <span class="dt">caption =</span> <span class="st">"Parameters diagnostics"</span>
)</code></pre></div>
<table class="table">
<caption>Parameters diagnostics</caption>
<thead><tr class="header">
<th align="left">model</th>
<th align="right">chain</th>
<th align="left">parameter</th>
<th align="right">mean</th>
<th align="right">se_mean</th>
<th align="right">sd</th>
<th align="right">X2.5.</th>
<th align="right">X25.</th>
<th align="right">X50.</th>
<th align="right">X75.</th>
<th align="right">X97.5.</th>
<th align="right">n_eff</th>
<th align="right">Rhat</th>
<th align="right">true</th>
<th align="right">rank</th>
<th align="right">n</th>
</tr></thead>
<tbody>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu11</td>
<td align="right">-10.10</td>
<td align="right">0.10</td>
<td align="right">1.14</td>
<td align="right">-12.15</td>
<td align="right">-10.96</td>
<td align="right">-10.04</td>
<td align="right">-9.26</td>
<td align="right">-8.13</td>
<td align="right">125.00</td>
<td align="right">1.01</td>
<td align="right">-11.47</td>
<td align="right">0.10</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu12</td>
<td align="right">-10.64</td>
<td align="right">0.01</td>
<td align="right">0.09</td>
<td align="right">-10.84</td>
<td align="right">-10.68</td>
<td align="right">-10.64</td>
<td align="right">-10.58</td>
<td align="right">-10.47</td>
<td align="right">125.00</td>
<td align="right">0.99</td>
<td align="right">-9.91</td>
<td align="right">1.00</td>
<td align="right">1</td>
</tr>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu21</td>
<td align="right">-0.03</td>
<td align="right">0.10</td>
<td align="right">1.01</td>
<td align="right">-1.87</td>
<td align="right">-0.77</td>
<td align="right">-0.12</td>
<td align="right">0.70</td>
<td align="right">2.00</td>
<td align="right">92.91</td>
<td align="right">0.99</td>
<td align="right">-0.46</td>
<td align="right">0.34</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu22</td>
<td align="right">9.74</td>
<td align="right">0.01</td>
<td align="right">0.13</td>
<td align="right">9.47</td>
<td align="right">9.67</td>
<td align="right">9.75</td>
<td align="right">9.81</td>
<td align="right">9.96</td>
<td align="right">125.00</td>
<td align="right">0.99</td>
<td align="right">0.75</td>
<td align="right">0.00</td>
<td align="right">1</td>
</tr>
<tr class="odd">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu31</td>
<td align="right">9.91</td>
<td align="right">0.12</td>
<td align="right">0.95</td>
<td align="right">8.23</td>
<td align="right">9.28</td>
<td align="right">9.72</td>
<td align="right">10.58</td>
<td align="right">11.98</td>
<td align="right">62.18</td>
<td align="right">1.00</td>
<td align="right">10.48</td>
<td align="right">0.74</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">Dummy Model</td>
<td align="right">1</td>
<td align="left">mu32</td>
<td align="right">0.57</td>
<td align="right">0.01</td>
<td align="right">0.15</td>
<td align="right">0.28</td>
<td align="right">0.47</td>
<td align="right">0.58</td>
<td align="right">0.68</td>
<td align="right">0.85</td>
<td align="right">125.00</td>
<td align="right">1.00</td>
<td align="right">9.60</td>
<td align="right">1.00</td>
<td align="right">1</td>
</tr>
</tbody>
</table>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class="kw">plot</span>(
  val<span class="op">$</span>parameters<span class="op">$</span>Rhat, 
  <span class="dt">type =</span> <span class="st">"h"</span>, <span class="dt">ylim =</span> <span class="kw">c</span>(<span class="fl">0.9</span>, <span class="fl">1.15</span>),
  <span class="dt">ylab =</span> <span class="kw">expression</span>(<span class="kw">hat</span>(R)),
  <span class="dt">xlab =</span> <span class="st">"Validation run (n)"</span>
)
<span class="kw">abline</span>(<span class="dt">h =</span> <span class="fl">1.1</span>)</code></pre></div>
<p><img src="introduction_files/figure-html/unnamed-chunk-9-1.png" width="700"></p>
</div>
<div id="simulate-fit" class="section level2">
<h2 class="hasAnchor">
<a href="#simulate-fit" class="anchor"></a>3. Simulate &amp; Fit</h2>
<p>Once specified, the user may:</p>
<ul>
<li>Call <code><a href="../reference/sim.html">sim()</a></code> to simulate <span class="math inline">\(N\)</span> datasets containing the observation variables <span class="math inline">\({\mathbf{y}}\)</span> and the hidden paths <span class="math inline">\(z\)</span>.</li>
<li>Call <code><a href="../reference/compile.html">compile()</a></code> to compile the underlying Stan model. The returned object may be later used for either sampling (MCMC) or optimization (MAP).
<ul>
<li>Call <code>sampling()</code> to draw posterior samples from a compiled model effectively allowing for full Bayesian inference via MCMC.</li>
<li>Call <code><a href="../reference/optimizing.html">optimizing()</a></code> to maximize the joint posterior density and obtain a point estimate.</li>
</ul>
</li>
<li>Call <code><a href="../reference/fit.html">fit()</a></code> to compile and sample in one step.</li>
</ul>
<p>After fitting, the underlying Stan code may be inspected by calling <code>browse_model</code>, which opens the Stan file in the editor. None that this function is only for reading purposes: changes in the Stan code will not affect the specification object.</p>
<div id="fully-bayesian-inferencia-via-mcmc" class="section level3">
<h3 class="hasAnchor">
<a href="#fully-bayesian-inferencia-via-mcmc" class="anchor"></a>3.1. Fully Bayesian Inferencia via MCMC</h3>
<p>Because <code>sampling()</code> and <code><a href="../reference/fit.html">fit()</a></code> return a <code>stanfit</code> S4 object, the posterior samples are compatible out of box with the <code>rstan</code> ecosystem [TRY AND LIST SOME SUCH AS bayesplot].</p>
</div>
<div id="maximum-a-posteriori-estimates-via-posterior-optimization" class="section level3">
<h3 class="hasAnchor">
<a href="#maximum-a-posteriori-estimates-via-posterior-optimization" class="anchor"></a>3.2. Maximum a Posteriori estimates via posterior optimization</h3>
<p><code>optimizing</code> is able to run several initialization of the optimization algorithm in parallel. The user may decide to <code>keep = "all"</code> runs, or just <code>keep = "best"</code>. If the former, <code>extract_grid</code> returns a summary of the different runs and <code>extract_best</code> retrieves the one with highest posterior log likelihood. If no run achieves convergence, a warning message is prompt to the console.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">1</span>,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="op">-</span><span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>   <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>  <span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">name =</span> <span class="st">"Univariate Gaussian"</span>
)

<span class="kw">set.seed</span>(<span class="dv">9000</span>)
y &lt;-<span class="st"> </span><span class="kw">as.matrix</span>(
  <span class="kw">c</span>(<span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="op">-</span><span class="dv">10</span>, <span class="dv">1</span>), <span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="dv">0</span>, <span class="dv">1</span>), <span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="dv">10</span>, <span class="dv">1</span>))
)
y &lt;-<span class="st"> </span><span class="kw">sample</span>(y)

myModel &lt;-<span class="st"> </span><span class="kw"><a href="../reference/compile.html">compile</a></span>(mySpec)
myAll   &lt;-<span class="st"> </span><span class="kw"><a href="../reference/optimizing.html">optimizing</a></span>(mySpec, myModel, <span class="dt">y =</span> y, <span class="dt">nRun =</span> <span class="dv">10</span>, <span class="dt">keep =</span> <span class="st">"all"</span>, <span class="dt">nCores =</span> <span class="dv">4</span>)
myBest &lt;-<span class="st"> </span><span class="kw"><a href="../reference/extract_best.html">extract_best</a></span>(myAll)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">print</span>(
  <span class="kw">round</span>(<span class="kw"><a href="../reference/extract_grid.html">extract_grid</a></span>(myAll, <span class="dt">pars =</span> <span class="st">"mu"</span>), <span class="dv">2</span>)
)
<span class="co">#&gt;        n       seed logPosterior returnCode         mu11  mu21  mu31</span>
<span class="co">#&gt;  [1,]  8 1111034741      -820.72          0 0 0 0 -10.06  9.98  0.04</span>
<span class="co">#&gt;  [2,]  2  768003877      -820.72          0 0 0 0 -10.06  9.98  0.04</span>
<span class="co">#&gt;  [3,]  6  820771557      -820.72          0 0 0 0 -10.06  9.98  0.04</span>
<span class="co">#&gt;  [4,] 10 2115105852      -820.72          0 0 0 0 -10.06  9.98  0.04</span>
<span class="co">#&gt;  [5,]  5 1539977502     -1018.84          0 0 0 0  -0.16  9.98 -9.87</span>
<span class="co">#&gt;  [6,]  9  355282361     -1018.84          0 0 0 0  -0.16  9.98 -9.86</span>
<span class="co">#&gt;  [7,]  7 2066428059     -1020.17          0 0 0 0   9.88 -9.96  0.04</span>
<span class="co">#&gt;  [8,]  3 1898572606     -1020.17          0 0 0 0   9.88 -9.96  0.04</span>
<span class="co">#&gt;  [9,]  1 1144725567     -1119.23          0 0 0 0   9.88 -0.06 -9.87</span>
<span class="co">#&gt; [10,]  4  433311989     -1119.23          0 0 0 0   9.88 -0.06 -9.87</span></code></pre></div>
<p><strong>NOTE</strong>: Passing both the specification and the compiled model is redundant. In future versions, user will not need to pass the specification.</p>
</div>
<div id="common-high-level-routines" class="section level3">
<h3 class="hasAnchor">
<a href="#common-high-level-routines" class="anchor"></a>3.3 Common high-level routines</h3>
<p>Most of these extractors rely on the lower level <code>extract_quantity</code> function. By default, it returns a named list where each element is an array representing a sampled quantity. The dimension of the array varies according to the number of iterations, the length of the series, and the structure of the quantity itself (for example, whether they are univariate or multivariate). Additionally, the function accepts wildcards <code>*</code> for easy extraction of a family of quantity.</p>
<p>This function has two convenience arguments:</p>
<ul>
<li>
<code>reduce</code> takes a function that is run with-in chain, effectively reducing the <span class="math inline">\(N\)</span> draws per chain into one or more quantity.</li>
<li>
<code>combine</code> takes a function that is run on the reduced quantity for each parameter, effectively combining the summary quantity from each parameter into a single data structure.</li>
</ul>
<p>Besides <code>extract_quantity</code>, which takes an argument <code>pars</code>, the software comes with extractors for the most typical quantities: <code>extract_alpha</code> (filtered probability), <code>extract_gamma</code> (smoothed probability), <code>extract_zstar</code> (MPM, also known a Viterbi), <code>extract_ypred</code> and <code>extract_zpred</code> for predictive quantities, and <code>extract_obs_parameters</code> and <code>extract_parameters</code> for observation model and all model parameters respectively.</p>
<p><strong>NOTE</strong>: Naming convention may be changed before submitting to CRAN. Consider for example <code>extract_alpha</code> versus <code>extract_filtered_probability</code>, or <code>extract_zstar</code> versus <code>extract_viterbi</code>.</p>
<p>Compare in these examples different configurations of <code>reduce</code> and <code>combine</code> :</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">1</span>,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="op">-</span><span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>   <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>  <span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">name =</span> <span class="st">"Univariate Gaussian"</span>
)

<span class="kw">set.seed</span>(<span class="dv">9000</span>)
y &lt;-<span class="st"> </span><span class="kw">as.matrix</span>(
  <span class="kw">c</span>(<span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="op">-</span><span class="dv">10</span>, <span class="dv">1</span>), <span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="dv">0</span>, <span class="dv">1</span>), <span class="kw">rnorm</span>(<span class="dv">100</span>, <span class="dv">10</span>, <span class="dv">1</span>))
)
y &lt;-<span class="st"> </span><span class="kw">sample</span>(y)

myFit &lt;-<span class="st"> </span><span class="kw"><a href="../reference/fit.html">fit</a></span>(mySpec, <span class="dt">y =</span> y, <span class="dt">iter =</span> <span class="dv">500</span>)

<span class="co"># 3 parameters, 1,000 iterations, 4 chains</span>
<span class="kw">print</span>(
  <span class="kw">str</span>(<span class="kw"><a href="../reference/extract_obs_parameters.html">extract_obs_parameters</a></span>(myFit))
)
<span class="co">#&gt; List of 3</span>
<span class="co">#&gt;  $ mu11: num [1:250, 1:4] 9.73 10.1 10.05 9.8 9.95 ...</span>
<span class="co">#&gt;  $ mu21: num [1:250, 1:4] -0.1933 -0.189 -0.1924 0.0366 -0.0616 ...</span>
<span class="co">#&gt;  $ mu31: num [1:250, 1:4] -9.93 -9.84 -9.82 -9.88 -9.77 ...</span>
<span class="co">#&gt; NULL</span>

<span class="co"># reduce within-chain draws to one quantity (median)</span>
<span class="kw">print</span>(
  <span class="kw"><a href="../reference/extract_obs_parameters.html">extract_obs_parameters</a></span>(myFit, <span class="dt">reduce =</span> median)
)
<span class="co">#&gt; $mu11</span>
<span class="co">#&gt; [1]  9.8797135 -0.1670912 -0.1545372  9.8801167</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu21</span>
<span class="co">#&gt; [1] -0.04011546  9.98843546  9.98835770 -9.96346881</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu31</span>
<span class="co">#&gt; [1] -9.8604877 -9.8669435 -9.8609835  0.0377377</span>

<span class="co"># reduce within-chain draws to two quantities (quantiles)</span>
<span class="kw">print</span>(
  <span class="kw"><a href="../reference/extract_obs_parameters.html">extract_obs_parameters</a></span>(
    myFit, <span class="dt">reduce =</span> <span class="kw"><a href="../reference/posterior_intervals.html">posterior_intervals</a></span>(<span class="kw">c</span>(<span class="fl">0.1</span>, <span class="fl">0.9</span>))
  )
)
<span class="co">#&gt; $mu11</span>
<span class="co">#&gt;           [,1]        [,2]        [,3]      [,4]</span>
<span class="co">#&gt; [1,]  9.749743 -0.28719452 -0.25653428  9.742132</span>
<span class="co">#&gt; [2,] 10.008315 -0.03632407 -0.04306187 10.014145</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu21</span>
<span class="co">#&gt;             [,1]      [,2]      [,3]       [,4]</span>
<span class="co">#&gt; [1,] -0.19350995  9.856644  9.852452 -10.100755</span>
<span class="co">#&gt; [2,]  0.07658331 10.104012 10.112068  -9.824095</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu31</span>
<span class="co">#&gt;            [,1]      [,2]      [,3]        [,4]</span>
<span class="co">#&gt; [1,] -10.002818 -9.992977 -9.990071 -0.09496092</span>
<span class="co">#&gt; [2,]  -9.722103 -9.741113 -9.745545  0.17199537</span>

<span class="co"># combine quantities into a matrix</span>
<span class="kw">print</span>(
  <span class="kw"><a href="../reference/extract_obs_parameters.html">extract_obs_parameters</a></span>(
      myFit, 
      <span class="dt">reduce  =</span> median,
      <span class="dt">combine =</span> rbind
  )
)
<span class="co">#&gt;             [,1]       [,2]       [,3]       [,4]</span>
<span class="co">#&gt; mu11  9.87971351 -0.1670912 -0.1545372  9.8801167</span>
<span class="co">#&gt; mu21 -0.04011546  9.9884355  9.9883577 -9.9634688</span>
<span class="co">#&gt; mu31 -9.86048774 -9.8669435 -9.8609835  0.0377377</span></code></pre></div>
<p>Since classification is a frequent task in the data analysis workflow, we provide a family of <code>classify_*</code> functions that extract the estimated quantities and return the most likely state: <code>classify_alpha</code> and <code>classify_gamma</code> return the state with highest posterior probability at each time step, while <code>classify_zstar</code> returns the posterior mode of the jointly most likely state.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">print</span>(
  <span class="kw">str</span>(<span class="kw"><a href="../reference/classify_alpha.html">classify_alpha</a></span>(myFit))
)
<span class="co">#&gt;  int [1:300, 1:4] 3 2 2 1 1 2 2 1 1 2 ...</span>
<span class="co">#&gt; NULL</span>

<span class="kw">print</span>(
  <span class="kw">str</span>(<span class="kw"><a href="../reference/classify_gamma.html">classify_gamma</a></span>(myFit))
)
<span class="co">#&gt;  int [1:300, 1:4] 3 2 2 1 1 2 2 1 1 2 ...</span>
<span class="co">#&gt; NULL</span>

<span class="kw">print</span>(
  <span class="kw">str</span>(<span class="kw"><a href="../reference/classify_zstar.html">classify_zstar</a></span>(myFit))
)
<span class="co">#&gt;  num [1:4, 1:300] 3 3 3 2 2 1 1 3 2 1 ...</span>
<span class="co">#&gt; NULL</span></code></pre></div>
</div>
</div>
<div id="diagnose" class="section level2">
<h2 class="hasAnchor">
<a href="#diagnose" class="anchor"></a>4. Diagnose</h2>
<p>This step is not implemented as of the date of this writing. We plan to explore how to diagnose a HMM outside posterior predictive checks (think, for example, pseudo residuals).</p>
</div>
<div id="visualize" class="section level2">
<h2 class="hasAnchor">
<a href="#visualize" class="anchor"></a>5. Visualize</h2>
<p>We expect that visualization will play a very important role in our software once it arrives to a mature stage. At this point, we support a small amount of graphical routines:</p>
<ul>
<li>
<code>plot_series</code> focuses on the observation variables.</li>
<li>
<code>plot_state_probability</code> focuses on the estimated assignment probability.</li>
<li>
<code>plot_ppredictive</code> focuses on posterior predictive checks.</li>
</ul>
<p>Besides functions for fully-customizable individual plots such as <code>plot_ppredictive_density(y, yPred)</code> (note that the function takes as arguments the actual values and not a fitted object), the software includes a convenient way to automatically include one or more plots from the same family on a single layout. Each plot may toggle one or more “feature” such as colored lines, points, or backgrounds. Additionally, plots adjust automatically for multivariate observations, and variable names are automatically shown if the data matrix used to fitted the model has named columns. Colors are themeable via global options.</p>
<p>The following features are currently available:</p>
<ul>
<li><p>When plotting the observation variables: * <code>stateShade</code> to shade the background based on the state probability. * <code>yColoredLine</code> to include a line colored according to the state probability. * <code>yColoredDots</code> to include dots colored according to the state probability.</p></li>
<li><p>When plotting the state probabilities: * <code>stateShade</code> to shade the background based on the state probability. * <code>probabilityColoredLine</code> to include a line colored according to the state probability. * <code>probabilityColoredDots</code> to include dots colored according to the state probability. * <code>bottomColoredMarks</code> and <code>topColoredMarks</code> to include colored tick marks on the bottom or top axis. * <code>probabilityFan</code> to show the posterior interval of assignment probabilities.</p></li>
</ul>
<p>Again, the objects returned by <code>fit</code> and <code>sampling</code> are compatible with the rstan ecosystem, including visualization routines in the bayesplot package.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mySpec &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> <span class="dv">3</span>, <span class="dt">R =</span> <span class="dv">2</span>,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="op">-</span><span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>   <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span>  <span class="dv">10</span>, <span class="dt">sigma =</span> <span class="dv">1</span>), <span class="dt">sigma =</span> <span class="dv">1</span>),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">1.0</span>, <span class="fl">0.2</span>, <span class="fl">0.2</span>)) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.2</span>, <span class="fl">1.0</span>, <span class="fl">0.2</span>)) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.2</span>, <span class="fl">0.2</span>, <span class="fl">1.0</span>)),
  <span class="dt">name =</span> <span class="st">"Univariate Gaussian Dummy Model"</span>
)

mySim &lt;-<span class="st"> </span><span class="kw"><a href="../reference/sim.html">sim</a></span>(mySpec, <span class="dt">T =</span> <span class="dv">300</span>, <span class="dt">seed =</span> <span class="dv">9000</span>, <span class="dt">iter =</span> <span class="dv">500</span>, <span class="dt">chains =</span> <span class="dv">1</span>)
ySim  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/extract_ypred.html">extract_ypred</a></span>(mySim)[<span class="dv">1</span>, , ]
<span class="kw">colnames</span>(ySim) &lt;-<span class="st"> </span><span class="kw">c</span>(<span class="st">"Weight"</span>, <span class="st">"Height"</span>)

myFit &lt;-<span class="st"> </span><span class="kw"><a href="../reference/fit.html">fit</a></span>(mySpec, <span class="dt">y =</span> ySim, <span class="dt">seed =</span> <span class="dv">9000</span>, <span class="dt">iter =</span> <span class="dv">500</span>, <span class="dt">chains =</span> <span class="dv">1</span>)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/plot_series.html">plot_series</a></span>(myFit, <span class="dt">legend.cex =</span> <span class="fl">0.8</span>)

<span class="kw"><a href="../reference/plot_series.html">plot_series</a></span>(myFit, <span class="dt">xlab =</span> <span class="st">"Time steps"</span>, <span class="dt">features =</span> <span class="kw">c</span>(<span class="st">"yColoredLine"</span>))

<span class="kw"><a href="../reference/plot_series.html">plot_series</a></span>(
  myFit, <span class="dt">stateProbability =</span>  <span class="st">"smoothed"</span>, 
  <span class="dt">features =</span> <span class="kw">c</span>(<span class="st">"stateShade"</span>, <span class="st">"bottomColoredMarks"</span>)
)

<span class="kw"><a href="../reference/plot_state_probability.html">plot_state_probability</a></span>(myFit, <span class="dt">main =</span> <span class="st">"Title"</span>, <span class="dt">xlab =</span> <span class="st">"Time"</span>)

<span class="kw"><a href="../reference/plot_state_probability.html">plot_state_probability</a></span>(myFit, <span class="dt">features =</span> <span class="kw">c</span>(<span class="st">"stateShade"</span>), <span class="dt">main =</span> <span class="st">"Title"</span>, <span class="dt">xlab =</span> <span class="st">"Time"</span>)

<span class="co"># plot_state_probability(myFit, features = c("bottomColoredMarks"), main = "Title", xlab = "Time")</span>

<span class="kw"><a href="../reference/plot_state_probability.html">plot_state_probability</a></span>(myFit, <span class="dt">features =</span> <span class="kw">c</span>(<span class="st">"probabilityColoredDots"</span>), <span class="dt">main =</span> <span class="st">"Title"</span>, <span class="dt">xlab =</span> <span class="st">"Time"</span>)

<span class="kw"><a href="../reference/plot_state_probability.html">plot_state_probability</a></span>(myFit, <span class="dt">features =</span> <span class="kw">c</span>(<span class="st">"probabilityColoredLine"</span>), <span class="dt">main =</span> <span class="st">"Title"</span>, <span class="dt">xlab =</span> <span class="st">"Time"</span>)

<span class="co"># plot_state_probability(myFit, stateProbability = "filtered", features = c("bottomColoredMarks", "probabilityFan"), stateProbabilityInterval = c(0.05, 0.95), main = "Title", xlab = "Time")</span>

<span class="kw"><a href="../reference/plot_ppredictive.html">plot_ppredictive</a></span>(myFit, <span class="dt">type =</span> <span class="kw">c</span>(<span class="st">"density"</span>, <span class="st">"cumulative"</span>, <span class="st">"summary"</span>), <span class="dt">fun =</span> median)

<span class="kw"><a href="../reference/plot_ppredictive.html">plot_ppredictive</a></span>(myFit, <span class="dt">type =</span> <span class="kw">c</span>(<span class="st">"density"</span>, <span class="st">"boxplot"</span>), <span class="dt">fun =</span> median, <span class="dt">subset =</span> <span class="dv">1</span><span class="op">:</span><span class="dv">10</span>)

<span class="kw"><a href="../reference/plot_ppredictive.html">plot_ppredictive</a></span>(
  myFit, 
  <span class="dt">type =</span> <span class="kw">c</span>(<span class="st">"density"</span>, <span class="st">"boxplot"</span>, <span class="st">"scatter"</span>), 
  <span class="dt">fun =</span> median, <span class="dt">fun1 =</span> mean, <span class="dt">fun2 =</span> median, 
  <span class="dt">subset =</span> <span class="dv">1</span><span class="op">:</span><span class="dv">40</span>
)

<span class="kw"><a href="../reference/plot_ppredictive.html">plot_ppredictive</a></span>(myFit, <span class="dt">type =</span> <span class="kw">c</span>(<span class="st">"density"</span>, <span class="st">"cumulative"</span>, <span class="st">"hist"</span>), <span class="dt">fun =</span> median)

<span class="kw"><a href="../reference/plot_ppredictive.html">plot_ppredictive</a></span>(myFit, <span class="dt">type =</span> <span class="kw">c</span>(<span class="st">"density"</span>, <span class="st">"cumulative"</span>, <span class="st">"ks"</span>))</code></pre></div>
</div>
<div id="select" class="section level2">
<h2 class="hasAnchor">
<a href="#select" class="anchor"></a>6. Select</h2>
<p>Model comparison is an important step in the statistical workflow, and future versions of the software will provide information criterion and other quantities for model selection.</p>
<!--
# Mixture models

Mixture models may be seen as a simplified version of HMM, where the assignment probabilities do not have a Markovian structure. We provide a simple implementation, and future iterations of our software will include extracting the individual assignment probabilities, and support for visualization and other data analysis procedures specific to mixture models.




```r
# print_fit(myFit)
```

-->
</div>
</div>
<div id="acknowledgements" class="section level1">
<h1 class="hasAnchor">
<a href="#acknowledgements" class="anchor"></a>Acknowledgements</h1>
<p>We acknowledge <a href="https://summerofcode.withgoogle.com/projects/#4681157036212224">Google Summer of Code 2018</a> for funding.</p>
</div>
<div id="appendix" class="section level1">
<h1 class="hasAnchor">
<a href="#appendix" class="anchor"></a>Appendix</h1>

<p>We present 19 different possible configurations.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># OBSERVATION MODEL -------------------------------------------------------</span>

K =<span class="st"> </span><span class="dv">3</span>
R =<span class="st"> </span><span class="dv">2</span>

<span class="co"># Case 1. A different multivariate density for each state</span>
<span class="co">#   Input: K multivariate densities</span>
<span class="co">#   Behaviour: Nothing</span>

exCase1  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/MVGaussianCholeskyCor.html">MVGaussianCholeskyCor</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>),
      <span class="dt">L     =</span> <span class="kw">LKJCor</span>(<span class="dt">eta =</span> <span class="dv">2</span>)
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/MVGaussianCholeskyCor.html">MVGaussianCholeskyCor</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">10</span>),
      <span class="dt">L     =</span> <span class="kw">LKJCor</span>(<span class="dt">eta =</span> <span class="dv">3</span>)
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/MVGaussianCholeskyCor.html">MVGaussianCholeskyCor</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">100</span>),
      <span class="dt">L     =</span> <span class="kw">LKJCor</span>(<span class="dt">eta =</span> <span class="dv">4</span>)
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"A different multivariate density for each each state"</span>
)

<span class="co"># Case 2. Same multivariate density for every state</span>
<span class="co">#   Input: One multivariate density</span>
<span class="co">#   Behaviour: Repeat input K times</span>

exCase2  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/MVGaussianCholeskyCor.html">MVGaussianCholeskyCor</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">100</span>),
      <span class="dt">L     =</span> <span class="kw">LKJCor</span>(<span class="dt">eta =</span> <span class="dv">2</span>)
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same multivariate density for every state"</span>
)

<span class="co"># Case 3. Same univariate density for every state and every output variable</span>
<span class="co">#   Input: One univariate density</span>
<span class="co">#   Behaviour: Repeat input K %nested% R times</span>

exCase3  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same univariate density for every state and every output variable"</span>
)

<span class="co"># Case 4. Same R univariate densities for every state</span>
<span class="co">#   Input: R univariate densities</span>
<span class="co">#   Behaviour: Repeat input K times</span>

exCase4  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same R univariate densities for every state"</span>
)

<span class="co"># Case 5. Same univariate density for every output variable</span>
<span class="co">#   Input: K univariate densities</span>
<span class="co">#   Behaviour: Repeat input R times</span>

exCase5  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same univariate density for every output variable"</span>
)

<span class="co"># Case 6. Different univariate densities for every pair of state and output variable</span>
<span class="co">#   Input: K %nested% R univariate densities</span>
<span class="co">#   Behaviour: None</span>

exCase6  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Different univariate densities for every pair of state and output variable"</span>
)

<span class="co"># UNIVARIATE Observation model --------------------------------------------</span>

K =<span class="st"> </span><span class="dv">3</span>
R =<span class="st"> </span><span class="dv">1</span>

<span class="co"># Case 7. A different univariate density for each each state</span>
<span class="co">#   Input: K univariate densities</span>
<span class="co">#   Behaviour: Nothing</span>

exCase7  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"A different univariate density for each each state"</span>
)

<span class="co"># Case 8. Same univariate density for every state</span>
<span class="co">#   Input: One univariate density</span>
<span class="co">#   Behaviour: Repeat input K times</span>

exCase8  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial     =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same multivariate density for every state"</span>
)

<span class="co"># INITIAL MODEL -----------------------------------------------------------</span>
K =<span class="st"> </span><span class="dv">3</span>
R =<span class="st"> </span><span class="dv">2</span>

<span class="co"># Case 9. Same univariate density for every initial state</span>
<span class="co">#   Input: One univariate density</span>
<span class="co">#   Behaviour: Repeat input K times</span>
exCase9  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span>
    <span class="kw"><a href="../reference/Beta.html">Beta</a></span>(
      <span class="dt">alpha =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">1</span>),
      <span class="dt">beta  =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">1</span>, <span class="dv">10</span>)
    ),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"Same univariate density for every initial state"</span>
)

<span class="co"># Case 10. One multivariate density for the whole initial vector</span>
<span class="co">#   Input: One multivariate density</span>
<span class="co">#   Behaviour: Nothing</span>
exCase10  &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span>
    <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(
      <span class="dt">alpha =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>()
    ),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"One multivariate density for the whole initial vector"</span>
)

<span class="co"># Case 11. A different univariate density for each initial state</span>
<span class="co">#   Input: K univariate densities</span>
<span class="co">#   Behaviour: Nothing</span>
exCase11 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"A different univariate density for each initial state"</span>
)

<span class="co"># Case 12. A link</span>
<span class="co">#   Input: One link density</span>
<span class="co">#   Behaviour: Nothing</span>
exCase12 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span>
    <span class="kw"><a href="../reference/InitialSoftmax.html">InitialSoftmax</a></span>(
      <span class="dt">vBeta =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>()
    ),
  <span class="dt">transition  =</span> <span class="kw"><a href="../reference/ImproperUniform.html">ImproperUniform</a></span>(),
  <span class="dt">name =</span> <span class="st">"TV Initial distribution"</span>
)

<span class="co"># TRANSITION MODEL --------------------------------------------------------</span>
K =<span class="st"> </span><span class="dv">3</span>
R =<span class="st"> </span><span class="dv">2</span>

<span class="co"># Case 13. Same univariate density for every transition</span>
<span class="co">#   Input: One univariate density</span>
<span class="co">#   Behaviour: Repeat input KxK times</span>
exCase13 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>),
  <span class="dt">name =</span> <span class="st">"Same univariate density for every transition"</span>
)

<span class="co"># Case 14. Same multivariate density for every transition row</span>
<span class="co">#   Input: One multivariate density</span>
<span class="co">#   Behaviour: Repeat input K times</span>

exCase14 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(
      <span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.7</span>)
    ),
  <span class="dt">name =</span> <span class="st">"Same multivariate density for every transition row"</span>
)

<span class="co"># Case 15. A different univariate density for each element of the transition row</span>
<span class="co">#   Input: K univariate densities</span>
<span class="co">#   Behaviour: Repeat input K times</span>

exCase15 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.1</span>, <span class="dt">beta =</span> <span class="fl">0.1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.5</span>, <span class="dt">beta =</span> <span class="fl">0.5</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.9</span>, <span class="dt">beta =</span> <span class="fl">0.9</span>),
  <span class="dt">name =</span> <span class="st">"A different univariate density for each element of the transition row"</span>
)

<span class="co"># Case 16. A different multivariate density for each transition row</span>
<span class="co">#   Input: K multivariate densities</span>
<span class="co">#   Behaviour: nothing</span>

exCase16 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.1</span>, <span class="fl">0.1</span>, <span class="fl">0.1</span>)) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.9</span>, <span class="fl">0.9</span>, <span class="fl">0.9</span>)),
  <span class="dt">name =</span> <span class="st">"A different multivariate density for each transition row"</span>
)

<span class="co"># Case 17. Different univariate densities for each element of the transition matrix</span>
<span class="co">#   Input: KxK univariate densities</span>
<span class="co">#   Behaviour: None</span>

exCase17 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.1</span>, <span class="dt">beta =</span> <span class="fl">0.1</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.2</span>, <span class="dt">beta =</span> <span class="fl">0.2</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.3</span>, <span class="dt">beta =</span> <span class="fl">0.3</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.4</span>, <span class="dt">beta =</span> <span class="fl">0.4</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.5</span>, <span class="dt">beta =</span> <span class="fl">0.5</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.6</span>, <span class="dt">beta =</span> <span class="fl">0.6</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.7</span>, <span class="dt">beta =</span> <span class="fl">0.7</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.8</span>, <span class="dt">beta =</span> <span class="fl">0.8</span>) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Beta.html">Beta</a></span>(<span class="dt">alpha =</span> <span class="fl">0.9</span>, <span class="dt">beta =</span> <span class="fl">0.9</span>),
  <span class="dt">name =</span> <span class="st">"Different univariate densities for each element of the transition matrix"</span>
)

<span class="co"># Case 18. A link</span>
<span class="co">#   Input: One link density</span>
<span class="co">#   Behaviour: Nothing</span>

exCase18 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span> <span class="kw"><a href="../reference/Dirichlet.html">Dirichlet</a></span>(<span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>)),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/TransitionSoftmax.html">TransitionSoftmax</a></span>(
      <span class="dt">uBeta =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>)
    ),
  <span class="dt">name =</span> <span class="st">"Different univariate densities for each element of the transition matrix"</span>
)

<span class="co"># A FULLY COMPLEX MODEL ---------------------------------------------------</span>

<span class="co"># Case 19. A link in both the transition and initial distribution.</span>
<span class="co">#   Input: A link in both the transition and initial distribution.</span>
<span class="co">#   Behaviour: Nothing</span>
exCase19 &lt;-<span class="st"> </span><span class="kw"><a href="../reference/hmm.html">hmm</a></span>(
  <span class="dt">K =</span> K, <span class="dt">R =</span> R,
  <span class="dt">observation =</span>
    <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ) <span class="op">+</span>
<span class="st">    </span><span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(
      <span class="dt">mu    =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>),
      <span class="dt">sigma =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dv">0</span>, <span class="dv">10</span>, <span class="dt">bounds =</span> <span class="kw">list</span>(<span class="dv">0</span>, <span class="ot">NULL</span>))
    ),
  <span class="dt">initial =</span>
    <span class="kw"><a href="../reference/InitialSoftmax.html">InitialSoftmax</a></span>(
      <span class="dt">vBeta =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>)
    ),
  <span class="dt">transition  =</span>
    <span class="kw"><a href="../reference/TransitionSoftmax.html">TransitionSoftmax</a></span>(
      <span class="dt">uBeta =</span> <span class="kw"><a href="../reference/Gaussian.html">Gaussian</a></span>(<span class="dt">mu =</span> <span class="dv">0</span>, <span class="dt">sigma =</span> <span class="dv">1</span>)
    ),
  <span class="dt">name =</span> <span class="st">"Fully Complex Model"</span>
)</code></pre></div>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
        <div id="tocnav">
      <h2 class="hasAnchor">
<a href="#tocnav" class="anchor"></a>Contents</h2>
      <ul class="nav nav-pills nav-stacked">
<li><a href="#introduction">Introduction</a></li>
      <li>
<a href="#hidden-markov-models">Hidden Markov Models</a><ul class="nav nav-pills nav-stacked">
<li><a href="#observation-model">Observation model</a></li>
      <li><a href="#transition-model">Transition model</a></li>
      <li><a href="#initial-distribution-model">Initial distribution model</a></li>
      <li><a href="#inference">Inference</a></li>
      <li><a href="#parameter-estimation">Parameter estimation</a></li>
      <li><a href="#naming-convention">Naming convention</a></li>
      </ul>
</li>
      <li>
<a href="#using-bayeshmm">Using BayesHMM</a><ul class="nav nav-pills nav-stacked">
<li><a href="#specify">1. Specify</a></li>
      <li><a href="#validate">2. Validate</a></li>
      <li><a href="#simulate-fit">3. Simulate &amp; Fit</a></li>
      <li><a href="#diagnose">4. Diagnose</a></li>
      <li><a href="#visualize">5. Visualize</a></li>
      <li><a href="#select">6. Select</a></li>
      </ul>
</li>
      <li><a href="#acknowledgements">Acknowledgements</a></li>
      <li><a href="#appendix">Appendix</a></li>
      </ul>
</div>
      </div>

</div>


      <footer><div class="copyright">
  <p>Developed by Luis Damiano, Michael Weylandt, Brian Peterson.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="http://pkgdown.r-lib.org/">pkgdown</a>.</p>
</div>

      </footer>
</div>

  

  </body>
</html>
