2

我写了一些不起作用的小 pymc3 代码。 

import pymc3

def create_model_pymc(data):
    with pymc3.Model() as model:
        k = 3
        #discussion about conjugate prior of a dirichlet are not so common see http://andrewgelman.com/2009/04/29/conjugate_prior/
        #u = pymc3.Uniform("u", lower=.1, upper=data.max()+1., observed=data, shape=k)
        u = pymc3.Exponential('u', 1./10, observed=data, shape=k)
        p = pymc3.Dirichlet('p', a=u, shape=k)
        c = pymc3.Categorical('c', p=p)
    return model

alpha_posterior = np.array([10., 3., 4.])
model = create_model_pymc(alpha_posterior)

with model:
    step = pymc3.Metropolis(model.vars)
    trace = pymc3.sample(20000, step)
    #trace = pymc3.sample(20000)

我收到一个 theano 错误“IndexError: index out of bounds”错误消息，但我不知道为什么。我尝试了上面代码的不同变体，但它们都不起作用。

如果您对代码的含义感到好奇，这是为了在观察多项式计数后检查狄利克雷参数的后验：（10,3,4）对于具有“均匀”狄利克雷先验的三个类（alpha=[1 .,1.,1.]) 导致 (10.,3.,4.)

简单直接采样：

a=alpha_posterior
print a
nb_samples = 200000
c_ = np.ndarray(nb_samples)
for i in range(nb_samples):
    d = scipy.stats.dirichlet.rvs(a,1)[0]
    c_[i] = np.random.choice(3, 1, p=d)
ns = float(nb_samples)

print (c_==0).sum()/ns
print (c_==1).sum()/ns
print (c_==2).sum()/ns
print x
print x/float(x.sum())
4

1 回答 1

0

这段代码对我来说很好，没有错误。

import numpy as np
import pymc3

with pymc3.Model() as model:
    data = np.array([10., 3., 4.])
    k = 3
    u = pymc3.Exponential('u', 1./10, observed=data, shape=k)
    p = pymc3.Dirichlet('p', a=u, shape=k)
    c = pymc3.Categorical('c', p=p)
    step = pymc3.Metropolis()
    trace = pymc3.sample(20, step)

链中只有 20 个样本。

多进程采样（4 个作业中的 4 个链）

复合步骤

大都会：[c] 大都会：[p] 采样 4 条链：100%|█████████████████████████████████ ██████████████████████| 2080/2080 [00:02<00:00, 893.84draws/s]

于 2019-03-25T07:16:35.297 回答