为什么在这个逻辑回归示例中 Pymc3 ADVI 比 MCMC 差?
Why is Pymc3 ADVI worse than MCMC in this logistic regression example?
我知道 ADVI/MCMC 之间的数学差异,但我试图了解使用其中一个或另一个的实际含义。我正在 运行 以这种方式创建的数据上创建一个非常简单的逻辑回归示例:
import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
def logistic(x, b, noise=None):
L = x.T.dot(b)
if noise is not None:
L = L+noise
return 1/(1+np.exp(-L))
x1 = np.linspace(-10., 10, 10000)
x2 = np.linspace(0., 20, 10000)
bias = np.ones(len(x1))
X = np.vstack([x1,x2,bias]) # Add intercept
B = [-10., 2., 1.] # Sigmoid params for X + intercept
# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=0., size=len(x1)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
y = np.random.binomial(1., pnoisy)
而我 运行 ADVI 是这样的:
with pm.Model() as model:
# Define priors
intercept = pm.Normal('Intercept', 0, sd=10)
x1_coef = pm.Normal('x1', 0, sd=10)
x2_coef = pm.Normal('x2', 0, sd=10)
# Define likelihood
likelihood = pm.Bernoulli('y',
pm.math.sigmoid(intercept+x1_coef*X[0]+x2_coef*X[1]),
observed=y)
approx = pm.fit(90000, method='advi')
不幸的是,无论我增加多少采样,ADVI似乎都无法恢复我定义的原始贝塔[-10., 2., 1.],而MCMC工作正常(如下图)
感谢您的帮助!
这是个有趣的问题! PyMC3 中的默认 'advi'
是平均场变分推理,它不能很好地捕获相关性。原来你建立的模型有一个有趣的相关结构,可以看出:
import arviz as az
az.plot_pair(trace, figsize=(5, 5))
PyMC3 有一个内置的收敛检查器 - 运行 优化太长或太短会导致有趣的结果:
from pymc3.variational.callbacks import CheckParametersConvergence
with model:
fit = pm.fit(100_000, method='advi', callbacks=[CheckParametersConvergence()])
draws = fit.sample(2_000)
这对我来说在大约 60,000 次迭代后停止。现在我们可以检查相关性并看到,正如预期的那样,ADVI 适合轴对齐的高斯分布:
az.plot_pair(draws, figsize=(5, 5))
最后,我们可以比较 NUTS 和(平均场)ADVI 的拟合:
az.plot_forest([draws, trace])
请注意,ADVI 低估了方差,但与每个参数的平均值相当接近。此外,您可以设置 method='fullrank_advi'
以更好地捕获您看到的相关性。
(注意:arviz
即将成为 PyMC3 的绘图库)
我知道 ADVI/MCMC 之间的数学差异,但我试图了解使用其中一个或另一个的实际含义。我正在 运行 以这种方式创建的数据上创建一个非常简单的逻辑回归示例:
import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
def logistic(x, b, noise=None):
L = x.T.dot(b)
if noise is not None:
L = L+noise
return 1/(1+np.exp(-L))
x1 = np.linspace(-10., 10, 10000)
x2 = np.linspace(0., 20, 10000)
bias = np.ones(len(x1))
X = np.vstack([x1,x2,bias]) # Add intercept
B = [-10., 2., 1.] # Sigmoid params for X + intercept
# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=0., size=len(x1)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
y = np.random.binomial(1., pnoisy)
而我 运行 ADVI 是这样的:
with pm.Model() as model:
# Define priors
intercept = pm.Normal('Intercept', 0, sd=10)
x1_coef = pm.Normal('x1', 0, sd=10)
x2_coef = pm.Normal('x2', 0, sd=10)
# Define likelihood
likelihood = pm.Bernoulli('y',
pm.math.sigmoid(intercept+x1_coef*X[0]+x2_coef*X[1]),
observed=y)
approx = pm.fit(90000, method='advi')
不幸的是,无论我增加多少采样,ADVI似乎都无法恢复我定义的原始贝塔[-10., 2., 1.],而MCMC工作正常(如下图)
感谢您的帮助!
这是个有趣的问题! PyMC3 中的默认 'advi'
是平均场变分推理,它不能很好地捕获相关性。原来你建立的模型有一个有趣的相关结构,可以看出:
import arviz as az
az.plot_pair(trace, figsize=(5, 5))
PyMC3 有一个内置的收敛检查器 - 运行 优化太长或太短会导致有趣的结果:
from pymc3.variational.callbacks import CheckParametersConvergence
with model:
fit = pm.fit(100_000, method='advi', callbacks=[CheckParametersConvergence()])
draws = fit.sample(2_000)
这对我来说在大约 60,000 次迭代后停止。现在我们可以检查相关性并看到,正如预期的那样,ADVI 适合轴对齐的高斯分布:
az.plot_pair(draws, figsize=(5, 5))
最后,我们可以比较 NUTS 和(平均场)ADVI 的拟合:
az.plot_forest([draws, trace])
请注意,ADVI 低估了方差,但与每个参数的平均值相当接近。此外,您可以设置 method='fullrank_advi'
以更好地捕获您看到的相关性。
(注意:arviz
即将成为 PyMC3 的绘图库)