<commit_before>## Sampling Inference for Bayesian HSMMs ##
`pyhsmm` is a nascent Python library for approximate unsupervised sampling
inference in Bayesian Hidden Markov Models (HMMs) and explicit-duration Hidden
semi-Markov Models (HSMMs), focusing on the Bayesian Nonparametric extensions,
the HDP-HMM and HDP-HSMM, via the weak-limit approximation. 

<!--
In the Bayesian paradigm, inference refers to both what would in other contexts
be called "learning" (or "parameter fitting") as well as "inference": all the
latent variables in the model, including hidden states and transition/emission
parameters, are included in the posterior distribution. The goal of sampling
inference is to produce (approximate) samples from the posterior, and each
sample roughly represents an alternative HMM or HSMM to explain the data. Using
the Bayesian Nonparametric HDP-HMM and HDP-HSMM, the sampled models that come
out can be of different complexity: there may be good explanations using only 5
states as well as good explanations that use 15 states. The purpose of this
sampling code is to produce samples of those alternatives.
-->

<!--
(At some point in the not-too-distant future, `py-hsmm` may also support
approximate learning/fitting of HMMs and HSMMs through the Expectation
Maximization (EM) algorithm.)
-->

## Installing ##
You can clone this library with the usual command:

```
git clone git://github.com/mattjj/pyhsmm.git
```

There is an optional dependency on [Eigen C++ Template Library](http://eigen.tuxfamily.org/index.php?title=Main_Page) installed in the usual location of `/usr/local/include` (its default install location), which can make the inference run much faster in many cases (and only a bit faster in others). If you install Eigen, you can enable its use by calling `pyhsmm.use_eigen()` after importing `pyhsmm`.

## A Simple Demonstration ##
Here's how to draw from the HDP-HSMM posterior over HSMMs given a sequence of
observations. (The same example can be found in `examples/basic.py`.)

First, we'll generate a sequence of observations from a 4-state HSMM with
Poisson durations and 2D Gaussian emissions. Our observation sequence will be
500 samples long. The average state durations will be 10, 20, 30, and 40
samples. The details of the observation hyperparameters aren't important; we're
just sampling random 2D Gaussian distributions.

```python
# imports
import numpy as np
from matplotlib import pyplot as plt

import pyhsmm

# set parameters and hyperparameters
N = 4; T = 500; obs_dim = 2
durparams = (10,20,30,40)
obs_hypparams = {'mu_0':np.zeros(obs_dim),
                'lmbda_0':np.eye(obs_dim),
                'kappa_0':0.2,
                'nu_0':obs_dim+2}

# instantiate the distribution objects
true_obs_distns = [pyhsmm.observations.gaussian(**obs_hypparams) for state in xrange(N)]
true_dur_distns = [pyhsmm.durations.poisson(lmbda=param) for param in durparams]

# build the true HSMM
truemodel = pyhsmm.hsmm(truth_obs_distns,truth_dur_distns)
```

Next, we'll sample some data from the model and plot the model on top of the
data

```python
data, labels = truemodel.generate(T)
truemodel.plot()
plt.gcf().suptitle('True HSMM')
```

![Randomly-generated model and data](http://www.mit.edu/~mattjj/github/py-hsmm/truth.png)

Now that we've generated some data, we can try to infer the model that it came
from (not peeking into `truemodel`, of course; that would be cheating!).

First, we create a new HSMM object that doesn't know anything about the
original (though we'll give it the the same hyperparameters). We'll ask this
model to infer the number of states as well (since an HDP-HSMM is instantiated
by default), so we'll give it an `Nmax` parameter instead of telling it `N=4`:

```python
Nmax = 10
obs_distns = [pyhsmm.observations.gaussian(**obs_hypparams) for state in xrange(Nmax)]
dur_distns = [pyhsmm.durations.poisson() for state in xrange(Nmax)]

posteriormodel = pyhsmm.hsmm(obs_distns,dur_distns)
```

Then, we add the data we want to condition on:

```python
posteriormodel.add_data(data)
```

Now we run a resampling loop. For each iteration of the loop, all the latent
variables of the model will be resampled by Gibbs sampling steps, including the
transition matrix, the observation means and covariances, the duration
parameters, and the hidden state sequence. We'll plot the samples every several
iterations.

```python
plot_every = 50
for idx in progprint_xrange(101):
    if (idx % plot_every) == 0:
        posteriormodel.plot()
        plt.gcf().suptitle('inferred HSMM after %d iterations (arbitrary colors)' % idx)

    posteriormodel.resample()
```

![Sampled models](http://www.mit.edu/~mattjj/github/py-hsmm/posterior_animation.gif)

## Speed ##

HSMMs constitute a much more powerful model class than plain-old HMMs, and that
enhanced power comes with a computational price: each sampling iteration for an
HSMM is much slower than that of an HMM. But that price is often worthwhile if
you want to place priors on state durations or have the model learn rich
duration distributions. In addition, the increased cost of each iteration often
pays for itself, since HSMM samplers empirically seem to take fewer iterations
to converge than comparable HMM samplers.

Using my nothing-special i7-920 desktop machine and a NumPy/SciPy built against
Intel's MKL BLAS (which generally outperforms ATLAS for these vectorized
operations) along with `pyhsmm.use_eigen()`, here's how long the demo
iterations took:

```
In [1]: run demo
.........................  [  25/101,    0.13sec avg,    9.63sec ETA ]
.........................  [  50/101,    0.12sec avg,    5.98sec ETA ]
.........................  [  75/101,    0.12sec avg,    3.10sec ETA ]
.........................  [ 100/101,    0.12sec avg,    0.12sec ETA ]
.
   0.12sec avg,   12.06sec total
```

## Extending the Code ##
To add your own observation or duration distributions, implement the interfaces defined in `abstractions.py`.

## References ##
* <NAME> and <NAME>, "Bayesian Nonparametric Hidden Semi-Markov Models." arXiv:1203.1365v1

* <NAME> and <NAME>, "The Hierarchical Dirichlet Process Hidden Semi-Markov Model." 26th Conference on Uncertainty in Artificial Intelligence (UAI 2010), Avalon, California, July 2010.
<commit_msg>tweaked readme to add links to references<commit_after>## Sampling Inference for Bayesian HSMMs ##
`pyhsmm` is a nascent Python library for approximate unsupervised sampling
inference in Bayesian Hidden Markov Models (HMMs) and explicit-duration Hidden
semi-Markov Models (HSMMs), focusing on the Bayesian Nonparametric extensions,
the HDP-HMM and HDP-HSMM, via the weak-limit approximation. 

<!--
In the Bayesian paradigm, inference refers to both what would in other contexts
be called "learning" (or "parameter fitting") as well as "inference": all the
latent variables in the model, including hidden states and transition/emission
parameters, are included in the posterior distribution. The goal of sampling
inference is to produce (approximate) samples from the posterior, and each
sample roughly represents an alternative HMM or HSMM to explain the data. Using
the Bayesian Nonparametric HDP-HMM and HDP-HSMM, the sampled models that come
out can be of different complexity: there may be good explanations using only 5
states as well as good explanations that use 15 states. The purpose of this
sampling code is to produce samples of those alternatives.
-->

<!--
(At some point in the not-too-distant future, `py-hsmm` may also support
approximate learning/fitting of HMMs and HSMMs through the Expectation
Maximization (EM) algorithm.)
-->

## Installing ##
You can clone this library with the usual command:

```
git clone git://github.com/mattjj/pyhsmm.git
```

There is an optional dependency on [Eigen C++ Template Library](http://eigen.tuxfamily.org/index.php?title=Main_Page) installed in the usual location of `/usr/local/include` (its default install location), which can make the inference run much faster in many cases (and only a bit faster in others). If you install Eigen, you can enable its use by calling `pyhsmm.use_eigen()` after importing `pyhsmm`.

## A Simple Demonstration ##
Here's how to draw from the HDP-HSMM posterior over HSMMs given a sequence of
observations. (The same example can be found in `examples/basic.py`.)

First, we'll generate a sequence of observations from a 4-state HSMM with
Poisson durations and 2D Gaussian emissions. Our observation sequence will be
500 samples long. The average state durations will be 10, 20, 30, and 40
samples. The details of the observation hyperparameters aren't important; we're
just sampling random 2D Gaussian distributions.

```python
# imports
import numpy as np
from matplotlib import pyplot as plt

import pyhsmm

# set parameters and hyperparameters
N = 4; T = 500; obs_dim = 2
durparams = (10,20,30,40)
obs_hypparams = {'mu_0':np.zeros(obs_dim),
                'lmbda_0':np.eye(obs_dim),
                'kappa_0':0.2,
                'nu_0':obs_dim+2}

# instantiate the distribution objects
true_obs_distns = [pyhsmm.observations.gaussian(**obs_hypparams) for state in xrange(N)]
true_dur_distns = [pyhsmm.durations.poisson(lmbda=param) for param in durparams]

# build the true HSMM
truemodel = pyhsmm.hsmm(truth_obs_distns,truth_dur_distns)
```

Next, we'll sample some data from the model and plot the model on top of the
data

```python
data, labels = truemodel.generate(T)
truemodel.plot()
plt.gcf().suptitle('True HSMM')
```

![Randomly-generated model and data](http://www.mit.edu/~mattjj/github/py-hsmm/truth.png)

Now that we've generated some data, we can try to infer the model that it came
from (not peeking into `truemodel`, of course; that would be cheating!).

First, we create a new HSMM object that doesn't know anything about the
original (though we'll give it the the same hyperparameters). We'll ask this
model to infer the number of states as well (since an HDP-HSMM is instantiated
by default), so we'll give it an `Nmax` parameter instead of telling it `N=4`:

```python
Nmax = 10
obs_distns = [pyhsmm.observations.gaussian(**obs_hypparams) for state in xrange(Nmax)]
dur_distns = [pyhsmm.durations.poisson() for state in xrange(Nmax)]

posteriormodel = pyhsmm.hsmm(obs_distns,dur_distns)
```

Then, we add the data we want to condition on:

```python
posteriormodel.add_data(data)
```

Now we run a resampling loop. For each iteration of the loop, all the latent
variables of the model will be resampled by Gibbs sampling steps, including the
transition matrix, the observation means and covariances, the duration
parameters, and the hidden state sequence. We'll plot the samples every several
iterations.

```python
plot_every = 50
for idx in progprint_xrange(101):
    if (idx % plot_every) == 0:
        posteriormodel.plot()
        plt.gcf().suptitle('inferred HSMM after %d iterations (arbitrary colors)' % idx)

    posteriormodel.resample()
```

![Sampled models](http://www.mit.edu/~mattjj/github/py-hsmm/posterior_animation.gif)

## Speed ##

HSMMs constitute a much more powerful model class than plain-old HMMs, and that
enhanced power comes with a computational price: each sampling iteration for an
HSMM is much slower than that of an HMM. But that price is often worthwhile if
you want to place priors on state durations or have the model learn rich
duration distributions. In addition, the increased cost of each iteration often
pays for itself, since HSMM samplers empirically seem to take fewer iterations
to converge than comparable HMM samplers.

Using my nothing-special i7-920 desktop machine and a NumPy/SciPy built against
Intel's MKL BLAS (which generally outperforms ATLAS for these vectorized
operations) along with `pyhsmm.use_eigen()`, here's how long the demo
iterations took:

```
In [1]: run demo
.........................  [  25/101,    0.13sec avg,    9.63sec ETA ]
.........................  [  50/101,    0.12sec avg,    5.98sec ETA ]
.........................  [  75/101,    0.12sec avg,    3.10sec ETA ]
.........................  [ 100/101,    0.12sec avg,    0.12sec ETA ]
.
   0.12sec avg,   12.06sec total
```

## Extending the Code ##
To add your own observation or duration distributions, implement the interfaces defined in `abstractions.py`.

## References ##
* <NAME> and <NAME>, [Bayesian Nonparametric Hidden Semi-Markov Models](http://arxiv.org/abs/1203.1365). arXiv:1203.1365v1

* <NAME> and <NAME>, [The Hierarchical Dirichlet Process Hidden Semi-Markov Model](http://www.mit.edu/~mattjj/papers/uai2010.pdf). 26th Conference on Uncertainty in Artificial Intelligence (UAI 2010), Avalon, California, July 2010.

{'fraction_non_alphanumeric': 0.09716183097478062, 'fraction_numerical': 0.02023875405170369, 'mean_word_length': 4.1276854479124445, 'pattern_counts': {'":': 2, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '850028', 'n_tokens_mistral': 4091, 'n_tokens_neox': 3671, 'n_words': 1582}