0

I am trying use PYMC3 to implement an example where the data comes from a mixture of multinomials. The goal is to infer the underlying state_prob vector (see below). The code runs, but the Metropolis sampler gets stuck at the initial state_prior vector. Also, for some reason I have not been able to get NUTS to work.

import numpy as np
import pandas as pd

from pymc3 import Model, Multinomial, Dirichlet
import pymc3
import theano.tensor as tt
from theano import function, printing

N = 10
state_prior = np.array([.3, .3, .3])
state_prob = np.array([0.3, 0.1, 0.6]) #need to infer this
state_emission_tran = np.array([[0.3, 0.2, 0.5],
                               [0.1, 0.5, 0.4],
                               [0.0, 0.05, 0.95]])

state_data = np.random.multinomial(1, state_prob, size=N)
emission_prob_given_state = np.matmul(state_data, state_emission_tran)
def rand_mult(row_p):
    return np.random.multinomial(1, row_p)
emission_data = np.apply_along_axis(rand_mult, 1, emission_prob_given_state)
# done with creating data

with Model() as simple_dag:

    state = Dirichlet('state', state_prior*100, shape=3)
    emission_dist = [pymc3.Multinomial.dist(p=state_emission_tran[i], n=1, shape=3) for i in range(3)]
    emission_mix = pymc3.Mixture('emission_mix', w = state, comp_dists = emission_dist, observed=emission_data)

with simple_dag:
    step = pymc3.Metropolis(vars=[state])
    trace = pymc3.sample(10000, cores=2, chains=2, tune=500, step=step, progressbar=True)
Antony Joseph
  • 255
  • 2
  • 7

1 Answers1

0

Try this one:

import numpy as np
import pandas as pd

from pymc3 import Model, Multinomial, Dirichlet
import pymc3
import theano.tensor as tt
from theano import function, printing

N = 10
state_prior = np.array([.3, .3, .3])
state_prob = np.array([0.3, 0.1, 0.6]) #need to infer this
state_emission_tran = np.array([[0.3, 0.2, 0.5],
                               [0.1, 0.5, 0.4],
                               [0.0, 0.05, 0.95]])

state_data = np.random.multinomial(1, state_prob, size=N)
emission_prob_given_state = np.matmul(state_data, state_emission_tran)
def rand_mult(row_p):
    return np.random.multinomial(1, row_p)
emission_data = np.apply_along_axis(rand_mult, 1, emission_prob_given_state)
# done with creating data

with Model() as simple_dag:

    state = Dirichlet('state', state_prior*100, shape=3)
    emission_dist = [pymc3.Multinomial.dist(p=state_emission_tran[i], n=1, shape=3) for i in range(3)]
    emission_mix = pymc3.Mixture('emission_mix', w = state, comp_dists = emission_dist, observed=emission_data)

with simple_dag:
    trace = pymc3.sample(3000, tune=1000)

I am using pymc3 version 3.5 in Linux and it works fine.

David
  • 435
  • 3
  • 12