python之将 numpy 函数转换为 theano

我正在使用 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())来加快速度的任何想法？我已经使用了几年，所以不知道正确的方法。

请您参考如下方法:

不建议使用开关函数，因为它会破坏参数空间的良好几何形状，并使使用像 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 的讨论。