在pymc3中如何计算GP在一条记录线上的对数后验

Question

用例

假设我在X_0有一个观察y_0，我想用一个超参数的高斯过程来建模。假设我然后通过分层抽样后验点来确定超参数theta中的一个分布。在

现在，我想评估另一个观察的对数后验概率，比如y_1，在X_1处，在超参数分布上取平均值， E_theta [ log P(y_1 | y_0, X_0, X_1, theta) ] 理想情况下，我会从后面画出theta，然后计算log P(y_1 | y_0, X_0, X_1, theta)，然后取几何平均值。在

问题1

假设我有一个采样的结果，即一个跟踪，例如：

with pm.Model() as model:
    ...
    trace = pm.sample(1000)

如何计算这些样本（或其中的一个子集）上的另一个张量？我只能使用定义为模型一部分的pm.Deterministic找到部分解决方案

^{pr2}$

这感觉不太对劲。我觉得你应该能够在采样后对跟踪的子集进行评估。在

问题2

在pymc3中，pm.gp有没有一个方法部分可以创建表示log P(y_1 | y_0 X_0 X_1 theta)的张量，或者我必须自己创建它，它需要写出后验均值和协方差（在pm.gp中已经做过了），然后调用cholesky decomp等

Answer 1

我将用这个用例的一个例子来回答。 从本质上讲，把抽样和条件定义分开， 然后使用下面方法2的方法1。在

import numpy as np
import pymc3 as pm

# Data generation
X0 = np.linspace(0, 10, 100)[:,None]
y0 = X0**(0.5) + np.exp(-X0/5)*np.sin(10*X0)
y0 += 0.1*np.random.normal(size=y0.shape)
y0 = np.squeeze(y0)

# y1
X1 = np.linspace(0, 15, 200)[:,None]
y1 = X1**(0.5) + np.exp(-X1/6)*np.sin(8*X1)
y1 = np.squeeze(y1)

# 'Solve' the inference problem
with pm.Model() as model:
    l = pm.HalfNormal('l',5.)
    cov_func = pm.gp.cov.ExpQuad(1, ls=l)
    gp = pm.gp.Marginal(cov_func=cov_func)
    y0_ = gp.marginal_likelihood('y0',X0,y0,0.1)
    trace = pm.sample(100)

# Define the object P(y1 | X1 y0 X0)
with model:
    y1_ = gp.conditional('y1',X1,given={'X':X0,'y':y0,'noise':0.1})
    # Note the given=... is not strictly required as it's cached from above

###
# Method 1

logp = y1_.logp
logp_vals1 = []
for point in trace:
    point['y1'] = y1
    logp_vals1.append(logp(point))
    # note this is approximately 100x faster than logp_vals1.append(y1_.logp(point))
    # because logp is a property with a lot of overhead

###
# Method 2

import theano
y1_shr = theano.shared(y1)
with model:
    logp = pm.Deterministic('logp', y1_.distribution.logp(y1_shr))

logp_val2 = [pm.distributions.draw_values([logp], point) for point in trace]

方法1在我的机器上似乎快了2-3倍。在