将numpy函数转换为theano

Question

我正在PyMC3计算一些我不会在这里讨论的内容，但是如果您有兴趣的话，可以从此链接中了解。

“ 2-lambdas”情况基本上是一个开关函数，需要将其编译为一个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

关于如何使用pure加快速度的任何想法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的话语。