pymc3中的多元线性回归
Multivariate linear regression in pymc3
在专门使用 emcee
多年后,我最近开始学习 pymc3
,但我 运行 遇到了一些概念性问题。
我正在练习 Hogg's Fitting a model to data 的第 7 章。这涉及到具有任意 2d 不确定性的直线的 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
所以我试图传递 MvNormal
x 和 y 的测量值 (mu
s) 及其相关的测量不确定性(y_err
和 x_err
) .但它似乎不喜欢 2d tau
参数。
有什么想法吗?这一定是可以的
谢谢
您可以尝试适配以下模型。是一个"regular"线性回归。但是 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().
效果好得多
您可能还想查看此 repository by David Hogg. And the book Statistical Rethinking,其中 McElreath 讨论了进行线性回归的问题,包括输入和输出变量中的错误。
在专门使用 emcee
多年后,我最近开始学习 pymc3
,但我 运行 遇到了一些概念性问题。
我正在练习 Hogg's Fitting a model to data 的第 7 章。这涉及到具有任意 2d 不确定性的直线的 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
所以我试图传递 MvNormal
x 和 y 的测量值 (mu
s) 及其相关的测量不确定性(y_err
和 x_err
) .但它似乎不喜欢 2d tau
参数。
有什么想法吗?这一定是可以的
谢谢
您可以尝试适配以下模型。是一个"regular"线性回归。但是 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().
效果好得多您可能还想查看此 repository by David Hogg. And the book Statistical Rethinking,其中 McElreath 讨论了进行线性回归的问题,包括输入和输出变量中的错误。