在PyMC中获取确定性变量的统计信息

Question

假设我有一组随机的(X,Y)点：

import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
import scipy

x = np.array(range(0,50))
y = np.random.uniform(low=0.0, high=40.0, size=200)
y = map((lambda a: a[0] + a[1]), zip(x,y))

plt.scatter(x,y)

                   

然后我用简单的线性回归来拟合这些点：

std = 20.
tau=1/(std**2)
alpha = pm.Normal('alpha', mu=0, tau=tau)
beta  = pm.Normal('beta', mu=0, tau=tau)
sigma = pm.Uniform('sigma', lower=0, upper=20)

y_est = alpha + beta * x

likelihood = pm.Normal('y', mu=y_est, tau=1/(sigma**2), observed=True, value=y)

model = pm.Model([likelihood, alpha, beta, sigma, y_est])
mcmc  = pm.MCMC(model)
mcmc.sample(40000, 15000)

我该如何得到 y_est[0]、y_est[1]、y_est[2] 等的分布或统计数据呢？（注意，这些变量对应的是每个输入 x 值的 估计 的 y 值。）

Answer 1

在PyMC 2中，如果你想查看一个确定性变量的轨迹，你需要把这个确定性变量放在一个Lambda对象里（或者用@deterministic来装饰一个函数）。在你的例子中，可以这样做：

y_est = Lambda('y_est', lambda a=alpha, b=beta: a + b * x)

然后你就可以像对待随机变量一样，调用summary方法或者绘制这个节点的图。

顺便说一下，你不需要自己创建一个Model对象，因为MCMC已经为你处理好了。你只需要：

mcmc = pm.MCMC([likelihood, alpha, beta, sigma, y_est])

甚至可以更简洁地写：

mcmc = pm.MCMC(vars())

Answer 2

根据@Chris的建议，下面的代码可以正常工作：

x     = pm.Uniform('x', lower=xmin, upper=xmax)
alpha = pm.Normal('alpha', mu=0, tau=tau)
beta  = pm.Normal('beta', mu=0, tau=tau)
sigma = pm.Uniform('sigma', lower=0, upper=20)

# The deterministic:
y_gen = pm.Lambda('y_gen', lambda a=alpha, x=x, b=beta: a + b * x)

然后可以按照以下方式从中提取样本：

mcmc = pm.MCMC([x, y_gen])
mcmc.sample(n_total_samples, n_burn_in)

x_trace = mcmc.trace('x')[:]
y_trace = mcmc.trace('y_gen')[:]