From 7978d13de81684a234733aa2791dec91aac02950 Mon Sep 17 00:00:00 2001 From: Caleb Weinreb Date: Tue, 5 Sep 2023 09:41:14 -0400 Subject: [PATCH] added test for parallel hmm sampler --- dynamax/hidden_markov_model/inference_test.py | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/dynamax/hidden_markov_model/inference_test.py b/dynamax/hidden_markov_model/inference_test.py index f72e8a58..9a4babfe 100644 --- a/dynamax/hidden_markov_model/inference_test.py +++ b/dynamax/hidden_markov_model/inference_test.py @@ -2,6 +2,7 @@ import itertools as it import jax.numpy as jnp import jax.random as jr +from jax import vmap import dynamax.hidden_markov_model.inference as core import dynamax.hidden_markov_model.parallel_inference as parallel @@ -285,3 +286,32 @@ def test_parallel_smoother(key=0, num_timesteps=100, num_states=3): posterior = core.hmm_smoother(initial_probs, transition_matrix, log_likelihoods) posterior2 = parallel.hmm_smoother(initial_probs, transition_matrix, log_likelihoods) assert jnp.allclose(posterior.smoothed_probs, posterior2.smoothed_probs, atol=1e-1) + + +def test_parallel_posterior_sample( + key=0, num_timesteps=5, num_states=2, eps=1e-3, + num_samples=1000000, num_iterations=5 +): + if isinstance(key, int): + key = jr.PRNGKey(key) + + max_unique_size = 1 << num_timesteps + + def iterate_test(key_iter): + keys_iter = jr.split(key_iter, num_samples) + args = random_hmm_args(key_iter, num_timesteps, num_states) + + # Sample sequences from posterior + state_seqs = vmap(parallel.hmm_posterior_sample, (0, None, None, None), (0, 0))(keys_iter, *args)[1] + unique_seqs, counts = jnp.unique(state_seqs, axis=0, size=max_unique_size, return_counts=True) + blj_sample = counts / counts.sum() + + # Compute joint probabilities + blj = jnp.exp(big_log_joint(*args)) + blj = jnp.ravel(blj / blj.sum()) + + # Compare the joint distributions + return jnp.allclose(blj_sample, blj, rtol=0, atol=eps) + + keys = jr.split(key, num_iterations) + assert iterate_test(keys[0]) \ No newline at end of file