probcomp
HMM forward filtering and backward sampling
VKM: "HMM forward filtering + backward sampling, on a real-world example" https://github.com/probcomp/metaprob-clojure/milestones/5
Regarding forward filtering backward sampling -- it is two stage algorithm. Here are some tips:
Say we have an HMM with T time steps and N hidden states.
It first runs the forward algorithm in order to compute the marginal posterior distribution on the final hidden state p(x_T | y_1, .., y_T) [[ this takes the form of a discrete probability distribution on 1..N]]. For each t=1, .., T, we compute the "alpha" values iteratively \alpha_t(x_t) := p(x_t, y_1, ..., y_t) according to the pseudocode on Wikipedia, but using log-space (i.e. the log-sum-exp trick) to avoid underflow.
\alpha_1(x_1) = p(x_1) * p(y_1 | x_1) [[ initialization ]]
\alpha_t(x_t) = p(y_t | x_t) * \sum_{x_{t-1}=1}^N p(x_{t-1}, y_1, ..., y_{t-1}) * p(x_t | x_{t-1})
= p(y_t | x_t) * \sum_{x_{t-1}=1}^N \alpha_t(x_t) * p(x_t | x_{t-1}) for t=2, .., T
Then, given the alpha value for some t, we can compute p(x_t | y_1, .., y_t) by normalizing:
p(x_t | y_1, .., y_t) = \alpha_t(x_t) / \sum_{x_t=1}^N \alpha_t(x_t) for t = 1, .., T
Then, you start the "backward sampling pass" by sampling a concrete value for x_T from p(x_T | y_1, .., y_T), which was computed in the last step of the forward algorithm. Then, you iteratively sample a value for x_{t-1} from the discrete distribution p(x_{t-1} | x_t, y_{1:T}), which can be computed using the distribution p(x_t | y_{1:t}) that was computed in the forward pass, and the concrete value just sampled for x_t, for each t=T, .., 2.
p(x_{t-1} | x_t, y_{1:T}) = p(x_{t-1} | y_{1:t-1}) * p(x_t | x_{t-1}) / Z_t
where Z_t = \sum_{x_{t-1}=1}^N p(x_{t-1} | y_{1:t-1}) * p(x_t | x_{t-1})
The set of concrete values x_1, .., x_T generated during the backward sampling pass together constitute a joint sample from the distribution P(x_1, .., x_T | y_1, .., y_T).
Here is a Jupyter Notebook for HMM implemented in Venture; it might be useful: https://github.com/probcomp/Venturecxx/blob/20170623-schaechtle-new-optimization-routines/benchmarks/inference-benchmarks/hmm/demo.ipynb
Per @jar398, legacy code relating to this issue might be available in the python and converted directories at the legacy tag.
