首页 > 解决方案 > pymc3 多项式混合卡住了

问题描述

我正在尝试使用 PYMC3 来实现一个数据来自多项式混合的示例。目标是推断底层state_prob向量(见下文)。代码运行,但 Metropolis 采样器卡在初始state_prior向量上。此外,由于某种原因,我无法让 NUTS 工作。

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)

标签: pymcpymc3mcmc

解决方案


试试这个:

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)

pymc3在 Linux 中使用 3.5 版,它工作正常。


推荐阅读