# Fitting a Hidden Markov model with Pyro

Hi everybody, I am trying to fit an HMM (hidden markov model) with the pyro, a probabilistic programming library.

I generated a dataset to test my model with those rules

• hidden states : $$z_{t+1} = a z_t + b \epsilon_t$$ where $$a, b = 0.1, 0.1$$
• observable states: $$x_t = z_t + 0.1 \tilde{\epsilon}_t$$ where $$\epsilon_t. \tilde{\epsilon_t}$$ are gaussian noises with mean 0 and variance 1

Then, I want to infer the coeficients $$a$$ and $$b$$ with the variational inference tool of the Pyro library:

import numpy as np
import torch
import pyro
import pyro.distributions as dist
import torch.nn as nn

transition = lambda z: 0.1 * z - 1 +   0.1 * np.random.normal(0,1)
emission = lambda x:  x  + 0.1 * np.random.normal(0,1)

import pylab as plt
%matplotlib inline
def generate_one_sample(n=100):
z_tp = [0]
x_tp = [0]
for i in range(n):
z_tp.append(transition(z_tp[-1]))
x_tp.append(emission(z_tp[-1]))
return x_tp, z_tp
def generate_many_samples(n_samples=10, n_points=100):
return np.array([generate_one_sample(n_points)[1] for _ in range(n_samples)])

transition = lambda z: 0.1 * z - 1 +   0.1 * np.random.normal(0,1)
emission = lambda x:  x  + 0.1 * np.random.normal(0,1)

class Transition(nn.Module):
"""
Parameterizes the bernoulli observation likelihood p(x_t | z_t)
"""
def __init__(self):
super(Transition, self).__init__()
# initialize the three linear transformations used in the neural network
self.layer1= nn.Linear(1, 1)
self.layer2= nn.Linear(1, 1)
self.softplus = nn.Softplus()
def forward(self, z_t):
"""
Given the latent z at a particular time step t we return the vector of
probabilities ps that parameterizes the bernoulli distribution p(x_t|z_t)
"""
return self.layer1(z_t), self.softplus(self.layer2(z_t))

class DMM(nn.Module):
def __init__(self):
super(DMM, self).__init__()
self.transition = Transition().cuda()
self = self.cuda()
def model(self, data):
pyro.module("dmm", self)
n_sample = data.shape[0]
loc_z = 0 * torch.ones(n_sample).cuda()
scale_z =  torch.ones(n_sample).cuda()*0.1
for t in range(data.shape[1]):
if t != 0:
loc_z, scale_z = self.transition(z_t.reshape(-1, 1))
with pyro.poutine.scale(None, .1):
z_t = pyro.sample('z_{}'.format(t), pyro.distributions.Normal(loc_z, scale_z).independent(1))
loc_x = z_t; # self.emission(z_t.reshape(-1, 1))
x_t = pyro.sample('x_{}'.format(t), pyro.distributions.Normal(loc_x, .1).independent(1),
obs=data[:,t])

def guide(self, data):
pyro.module("dmm", self)
n_sample = data.shape[0]
loc_z = 0 * torch.ones(n_sample).cuda()
scale_z = torch.ones(n_sample).cuda()
for t in range(data.shape[1]):
z_t = pyro.sample('z_{}'.format(t), pyro.distributions.Normal(loc_z, scale_z).independent(1))
loc_z, scale_z = self.transition(z_t.reshape(-1, 1))
dmm = DMM().cuda()

adam_params = {"lr": 0.01, "betas": (0.95, 0.999)}
from pyro.infer import SVI, Trace_ELBO
svi = SVI(dmm.model, dmm.guide, optimizer, loss=Trace_ELBO())

import sys
torch.cuda.empty_cache
for i in range(500):
if i % 5000 == 0:
test_me_im_famous = torch.tensor(generate_many_samples(n_samples=100, n_points=7),dtype=torch.float32).cuda()
print('\n')
sys.stdout.write('epoch {},   {}    \r'.format(i, svi.step(test_me_im_famous)))
import pylab as plt
%matplotlib inline
x = np.arange(0,1, 0.01)
plt.plot(x,
dmm.transition(torch.arange(0,1,0.01).cuda().reshape(-1, 1))[0].cpu().flatten().detach().numpy(),
label='transition')
plt.plot(x, 0.1 * x - 1 , ls='--', label='real transition')
plt.legend(loc='best')


unfortunately after 500 iterations, I do not converge the correct coefficients. Has anybody already encountered such kind of problems ?