回归中估计后验预测

Question

假设我有一组随机的X和Y点：

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)

假设我用线性回归来把每个X值对应的Y值看作一个高斯分布（也就是正态分布），那么我该如何估计后验预测，也就是p(y|x)，对于每一个（可能的）X值呢？

有没有简单的方法可以用pymc或者scikit-learn来实现这个？

Answer 1

如果我理解你想要的没错的话，你可以使用PyMC的git版本（也就是PyMC3）和glm这个子模块来实现。比如说：

import numpy as np
import pymc as pm
import matplotlib.pyplot as plt 
from pymc import glm 

## Make some data
x = np.array(range(0,50))
y = np.random.uniform(low=0.0, high=40.0, size=50)
y = 2*x+y
## plt.scatter(x,y)

data = dict(x=x, y=y)
with pm.Model() as model:
    # specify glm and pass in data. The resulting linear model, its likelihood and 
    # and all its parameters are automatically added to our model.
    pm.glm.glm('y ~ x', data)
    step = pm.NUTS() # Instantiate MCMC sampling algorithm
    trace = pm.sample(2000, step)


##fig = pm.traceplot(trace, lines={'alpha': 1, 'beta': 2, 'sigma': .5});## traces
fig = plt.figure()
ax = fig.add_subplot(111)
plt.scatter(x, y, label='data')
glm.plot_posterior_predictive(trace, samples=50, eval=x,
                              label='posterior predictive regression lines')

这样你就能得到类似于这个的结果

你可能会对这些博客文章感兴趣： 1 和 2，我就是从这里获取的灵感。

补充说明 如果你想为每个x值获取y值，可以试试这个代码，我是通过研究glm的源代码找到的。

lm = lambda x, sample: sample['Intercept'] + sample['x'] * x ## linear model
samples=50 ## Choose to be the same as in plot call
trace_det = np.empty([samples, len(x)]) ## initialise
for i, rand_loc in enumerate(np.random.randint(0, len(trace), samples)):
    rand_sample = trace[rand_loc]
    trace_det[i] = lm(x, rand_sample)
y = trace_det.T
y[0]

抱歉如果这个方法不是最优雅的 - 希望你能理解其中的逻辑。