Gfyuki

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 ?

@robin-nicole consider asking this question on http://forum.pyro.ai where Pyro developers might be able to help out.