python - 将 numpy 函数转换为 theano

Question

我正在使用 PyMC3 来计算一些我不会在这里涉及的东西，但你可以从 this link 中得到这个想法。如果有兴趣。

'2-lambdas' case 基本上是一个 switch 函数，它需要编译成一个 Theano 函数以避免 dtype 错误，看起来像这样:

import theano
from theano.tensor import lscalar, dscalar, lvector, dvector, argsort

@theano.compile.ops.as_op(itypes=[lscalar, dscalar, dscalar], otypes=[dvector])
def lambda_2_distributions(tau, lambda_1, lambda_2):
        """
        Return values of `lambda_` for each observation based on the 
        transition value `tau`.
        """
        out = zeros(num_observations)
        out[: tau] = lambda_1  # lambda before tau is lambda1
        out[tau:] = lambda_2  # lambda after (and including) tau is lambda2
        return out

我试图将其概括为适用于“n-lambdas”，其中 taus.shape[0] = lambdas.shape[0] - 1，但我只能想出这个numpy 执行速度非常慢。

@theano.compile.ops.as_op(itypes=[lvector, dvector], otypes=[dvector])
def lambda_n_distributions(taus, lambdas):

    out = zeros(num_observations)
    np_tau_indices = argsort(taus).eval()
    num_taus = taus.shape[0]
    for t in range(num_taus):
        if t == 0:
            out[: taus[np_tau_indices[t]]] = lambdas[t]
        elif t == num_taus - 1:
            out[taus[np_tau_indices[t]]:] = lambdas[t + 1]
        else:
            out[taus[np_tau_indices[t]]: taus[np_tau_indices[t + 1]]] = lambdas[t]
    return out

关于如何使用纯 Theano(避免调用 .eval())来加快速度的任何想法？我已经使用了几年，所以不知道正确的方法。

Answer 1

不建议使用开关函数，因为它会破坏参数空间的良好几何形状，并使使用像 NUTS 这样的现代采样器进行采样变得困难。

相反，您可以尝试使用开关函数的连续松弛对其进行建模。这里的主要思想是将第一个切换点之前的速率建模为基线；并在每个切换点之后添加逻辑函数的预测:

def logistic(L, x0, k=500, t=np.linspace(0., 1., 1000)):
    return L/(1+tt.exp(-k*(t_-x0)))

with pm.Model() as m2:
    lambda0 = pm.Normal('lambda0', mu, sd=sd)
    lambdad = pm.Normal('lambdad', 0, sd=sd, shape=nbreak-1)
    trafo = Composed(pm.distributions.transforms.LogOdds(), Ordered())
    b = pm.Beta('b', 1., 1., shape=nbreak-1, transform=trafo,
                testval=[0.3, 0.5])
    theta_ = pm.Deterministic('theta', tt.exp(lambda0 +
                                          logistic(lambdad[0], b[0]) +
                                          logistic(lambdad[1], b[1])))
    obs = pm.Poisson('obs', theta_, observed=y)

    trace = pm.sample(1000, tune=1000)

我在这里也使用了一些技巧，例如，PyMC3 代码库中还没有的复合转换。您可以在此处查看完整代码:https://gist.github.com/junpenglao/f7098c8e0d6eadc61b3e1bc8525dd90d

如果您有更多问题，请发帖至https://discourse.pymc.io使用您的模型和(模拟)数据。我更经常地检查和回答 PyMC3 的讨论。