2

我正在尝试在 pymc3 中实现一个相当简单的模型。要点是我有一些数据是从一系列随机选择中生成的。选择可以被认为是多项式，并且该过程选择选择作为先前选择的函数。

类别的总体概率使用 Dirichlet 先验建模。

似然函数必须针对手头的数据进行定制。数据是过程输出的 0 和 1 列表。我已经在 pymc2 中成功制作了模型，你可以在这篇博文中找到它。这是一个为这个问题生成测试数据的python函数：

ps = [0.2,0.35,0.25,0.15,0.0498,1/5000]
def make(ps):
    out = []
    while len(out) < 5:
        n_spots = 5-len(out)
        sp = sum(ps[:n_spots+1])
        P = [x/sp for x in ps[:n_spots+1]]
        l = np.argwhere(np.random.multinomial(1,P)==1).ravel()[0]
        #if len(out) == 4:
        #    l = np.argwhere(np.random.multinomial(1,ps[:2])==1).ravel()[0]
        out.extend([1]*l)
        if (out and out[-1] == 1 and len(out) < 5) or l == 0:
            out.append(0)
        #print n_spots, l, len(out)
    assert len(out) == 5
    return out

当我学习/迁移到 pymc3 时，我正在尝试将观察到的数据输入到自定义似然函数中，并且在此过程中遇到了几个问题。这可能是因为这是我第一次使用 Theano，但我希望有人能给一些建议。

这是我的代码（使用上面的 make 函数）：

import numpy as np
import pymc3 as pm
from scipy import optimize
import theano.tensor as T
from theano.compile.ops import as_op
from collections import Counter

# This function gets the attributes of the data that are relevant for calculating the likelihood 
def scan(value):
    groups = []
    prev = False
    s = 0
    for i in xrange(5):
        if value[i] == 0:
            if prev:
                groups.append((s,5-(i-s)))
                prev = False
                s = 0
            else:
                groups.append((0,5-i))
        else:
            prev = True
            s += 1
    if prev:
        groups.append((s,4-(i-s)))
    return groups

# The likelihood calculation for a single data point
def like1(v,p):
    l = 1
    groups = scan(v)
    for n, s in groups:
        l *= p[n]/p[:s+1].sum()
    return T.log(l)

# my custom likelihood class
class CustomDist(pm.distributions.Discrete):
    def __init__(self, ps, data, *args, **kwargs):
        super(CustomDist, self).__init__(*args, **kwargs)
        self.ps = ps
        self.data = data

    def logp(self,v):
        all_l = 0
        for v, k in self.data.items():
            l = like1(v,self.ps)
            all_l += l*k
        return all_l

# model creation
model = pm.Model()

with model:
    probs = pm.Dirichlet('probs',a=np.array([0.5]*6),shape=6,testval=np.array([1/6.0]*6))
    output = CustomDist("rolls",ps=probs,data=data,observed=True)

我可以在大约一分钟左右找到 MAP（我的机器是 Windows 7，i7-4790 @3.6GHz）。MAP 与输入概率向量匹配良好，这至少意味着模型正确链接。

但是，当我尝试进行跟踪时，我的内存使用量猛增（高达几个 gig），实际上我还没有足够的耐心让模型完成编译。在跟踪之前，我已经等了 10 分钟 + NUTS 或 HMC 编译。但是，大都会步进工作得很好（并且比 pymc2 快得多）。

我是不是太希望 Theano 能够很好地处理非 theano 数据的 for 循环？有没有更好的方法来编写这段代码，以便 Theano 可以很好地使用它，或者我是否因为我的数据是自定义 python 类型而无法使用数组/矩阵操作进行分析而受到限制？

提前感谢您的建议和反馈。请让我知道可能需要澄清的内容！

4

0 回答 0