pymc3中的多元线性回归

Question

pymc3在专门使用了多年之后，我最近开始学习emcee，并且遇到了一些概念问题。

我正在练习霍格的《将模型拟合到数据》的第 7 章。这涉及到具有任意二维不确定性的直线的 mcmc 拟合。我在 中很容易地完成了这一点emcee，但pymc给我带来了一些问题。

它本质上归结为使用多元高斯似然。

这是我到目前为止所拥有的。

from pymc3 import  *

import numpy as np
import matplotlib.pyplot as plt

size = 200
true_intercept = 1
true_slope = 2

true_x = np.linspace(0, 1, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * true_x
# add noise

# here the errors are all the same but the real world they are usually not!
std_y, std_x = 0.1, 0.1 
y = true_regression_line + np.random.normal(scale=std_y, size=size)
x = true_x + np.random.normal(scale=std_x, size=size)

y_err = np.ones_like(y) * std_y
x_err = np.ones_like(x) * std_x

data = dict(x=x, y=y)

with Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
    # Define priors
    intercept = Normal('Intercept', 0, sd=20)
    gradient = Normal('gradient', 0, sd=20)


    # Define likelihood
    likelihood = MvNormal('y', mu=intercept + gradient * x,
                        tau=1./(np.stack((y_err, x_err))**2.), observed=y)

    # start the mcmc!
    start = find_MAP() # Find starting value by optimization
    step = NUTS(scaling=start) # Instantiate MCMC sampling algorithm
    trace = sample(2000, step, start=start, progressbar=False) # draw 2000 posterior samples using NUTS sampling

这会引发错误：LinAlgError: Last 2 dimensions of the array must be square

因此，我尝试传递MvNormalx 和 y ( s) 的测量值mu及其相关的测量不确定性 (y_err和x_err)。但它似乎不喜欢 2dtau参数。

有任何想法吗？这一定是可能的

谢谢

Answer 1

您可以尝试调整以下模型。是一个“常规”线性回归。但x和y已被高斯分布取代。在这里，我不仅假设输入和输出变量的测量值，而且还假设其误差的可靠估计（例如由测量设备提供）。如果您不相信这些错误值，您可以尝试从数据中估计它们。

with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)
    obs_x = pm.Normal('obs_x', mu=x, sd=x_err, shape=len(x))
    obs_y = pm.Normal('obs_y', mu=y, sd=y_err, shape=len(y))

    likelihood = pm.Normal('y', mu=intercept + gradient * obs_x,
                    sd=epsilon, observed=obs_y)

    trace = pm.sample(2000)

如果您根据数据估计误差，则可以合理地假设它们可能是相关的，因此，您可以使用多元高斯，而不是使用两个单独的高斯。在这种情况下，您最终会得到如下所示的模型：

df_data = pd.DataFrame(data)
cov = df_data.cov()

with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)

    obs_xy = pm.MvNormal('obs_xy', mu=df_data, tau=pm.matrix_inverse(cov), shape=df_data.shape)

    yl = pm.Normal('yl', mu=intercept + gradient * obs_xy[:,0],
                    sd=epsilon, observed=obs_xy[:,1])

mu, sds, elbo = pm.variational.advi(n=20000)
step =  pm.NUTS(scaling=model.dict_to_array(sds), is_cov=True)
trace = pm.sample(1000, step=step, start=mu)

请注意，在之前的模型中，协方差矩阵是根据数据计算的。如果你打算这样做，那么我认为最好使用第一个模型，但如果你要估计协方差矩阵，那么第二个模型可能是一个明智的方法。

对于第二个模型，我使用 ADVI 来初始化它。ADVI 是初始化模型的好方法，通常它比 find_MAP() 效果好得多。

您可能还想查看 David Hogg 的这个存储库。McElreath 在《统计反思》一书中讨论了线性回归的问题，包括输入和输出变量的误差。