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)