0
$\begingroup$

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()
from pyro.optim import Adam, ClippedAdam

adam_params = {"lr": 0.01, "betas": (0.95, 0.999)}
optimizer = Adam(adam_params)
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 ?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.