您可以尝试调整以下模型。是一个“常规”线性回归。但x
and 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() 效果好得多。
您可能还想检查一下存储库 https://github.com/davidwhogg/DataAnalysisRecipes作者:大卫·霍格。还有书统计反思 http://xcelab.net/rm/statistical-rethinking/其中 McElreath 讨论了线性回归的问题,包括输入和输出变量的误差。