NUTS 和 Metropolis 层次概率模型的收敛问题
Convergence issues on hierarchical probit model with NUTS and Metropolis
我正在尝试通过添加概率转换并使结果服从伯努利分布，将层次模型从 Gelman 和 Hill reproduced in PyMC3 here 扩展到二进制结果数据。现在我正在使用玩具数据，所以我知道真正的价值。 Alpha 应该是 .1，beta 应该是 .5。
模型 运行 在扩展前使用 NUTS 采样器没问题。一旦我添加它，估计值就会缓慢增加并持续增加，直到模型在 10 到 200 次迭代之间停滞不前。这是它一直到 120（相对较长 运行）时的图像。
在扩展之前，Metropolis 需要进行 200,000 次迭代才能很好地修复真实参数值，但它最终做到了。扩展后它停滞在 30k 和 50k 之间。与 NUTS 不同，当您尝试停止它时它会完全崩溃 post-stall，所以我没有图片。更早地停止它给出了一个大致超过零的 beta 估计值，但分布范围很广。
代码贴在下面。
不知道是采样问题，还是规格问题。有没有更好的方法来指定 Probit？关于其他采样器的任何提示可以尝试吗？我已经尽可能地剥离我的模型以进行测试，并在我添加概率扩展后将其缩小到破坏，但我不知道下一步该怎么做。
#Generate Data
n=100
#Determine how many observations per group
evts=np.random.randint(10,100,n)
#Determine which groups will be receive treatment
x_g=np.random.binomial(1,.3,n)
#pre-create distribution of betas for groups
mu = np.random.normal(.5,.2,n)
#preallocate space in a dataframe
i = np.zeros(evts.sum())
y_obs = pd.DataFrame({'y':i.copy(),
'x':i.copy(),
'grp':i.copy()},
index = range(evts.sum()))
#populate dataframe with simulated data
i=0
for grp in range(100):
#index of observations for a given group
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
#generate outcomes using
#different dgp depending on treatment
if x_g[grp] ==1:
#shortcut to make sure 1>p>0
p_i = max((.1 + mu[grp]),0.01)
p_i = min(p_i,1)
out = np.random.binomial(1,p_i,evts[grp])
else:
out = np.random.binomial(1,.1,evts[grp])
#Assign to dataframe
y_obs.loc[ind,'y'] = out
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs = y_obs.astype(int)
print('starting model')
with pm.Model() as test_model:
#hyperpriors
mu_a=pm.Normal('mu_a',mu=0, sd=100**2)
sig_a = pm.Uniform('sig_a',lower=0,upper=100)
mu_b=pm.Normal('mu_b',mu=0, sd=100**2)
sig_b = pm.Uniform('sig_b',lower=0,upper=100)
#priors
a = pm.Normal('a',mu=mu_a,sd = sig_a, shape=n)
b = pm.Normal('b',mu=mu_b,sd = sig_b, shape=n)
eps = pm.Uniform('eps',lower=0,upper=100)
est = a[y_obs.grp] + b[y_obs.grp] * y_obs.x
#I get correct estimates when I
#stop here using commented out line.
# y_hat = pm.Normal('y_hat',
# mu=est,
# sd=eps,
# observed = y_obs.y)
#Probit transformation:
y_hat = pm.Normal('y_hat',
mu=est,
sd=eps,
shape=y_obs.shape[0])
mu_y = tt.mean(y_hat)
eps_hat = tt.var(y_hat)
p_hat = 0.5 * (1 + tt.erf((y_hat-mu_y) / (eps_hat*tt.sqrt(2))))
y = pm.Bernoulli('y',p=p_hat, observed = y_obs.y)
with test_model:
#Either:
mu,sds,elbo = pm.variational.advi(n=100000)
step = pm.NUTS(scaling=test_model.dict_to_array(sds),
is_cov=True)
test_trace = pm.sample(200, step, start=mu)
#or
# step=pm.Metropolis()
# test_trace = pm.sample(50000)
pm.traceplot(test_trace)#[-5000::3])
注意：编辑以修复行中的错字：'step = pm.NUTS(scaling=test_model.dict_to_array(sds),`
编辑：我为最初 posted 模型的概率扩展制作了更好的模拟数据。 （原始数据生成是现在 ADVI 给出了更好的估计，所以它从正确的位置开始，但 NUTS 仍然很快停止（大约十次迭代）。Metropolis 直接失败：我做了第一轮 5000 次迭代，并得到尝试绘制轨迹时出错。
新数据生成：
n=100
evts=np.random.randint(10,100,n)
x_g=np.random.binomial(1,.3,n)
i = np.zeros(evts.sum())
mu = np.random.normal(.5,.2,n)
mu0 = np.random.normal(.1,.05,n)
y_obs = pd.DataFrame({'y':i.copy(),'x':i.copy(),'grp':i.copy()},index = range(evts.sum()))
i=0
for grp in range(100):
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
if x_g[grp] ==1:
est = mu0[grp] + mu[grp]
else:
est=mu0[grp]
p_hat = tt.nnet.sigmoid(est).eval()
y_obs.loc[ind,'y_hat'] = est
y_obs.loc[ind,'y'] = np.random.binomial(1,p_hat,len(ind))
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs['grp']=y_obs.grp.astype(np.int64)
当 pymc3 尝试绘制密度时大都市出错：
ValueError: v cannot be empty
也许我误解了你的意思，但这个模型不应该工作：
with pm.Model() as test_model:
#hyperpriors
mu_a = pm.Flat('mu_a')
sig_a = pm.HalfCauchy('sig_a', beta=2.5)
mu_b = pm.Flat('mu_b')
sig_b = pm.HalfCauchy('sig_b', beta=2.5)
#priors
a_raw = pm.Normal('a_raw', mu=0, sd=1, shape=n)
a = pm.Deterministic('a', mu_a + sig_a * a_raw)
b_raw = pm.Normal('b_raw', mu=0, sd=1, shape=n)
b = pm.Deterministic('b', mu_b + sig_b * b_raw)
est = a[y_obs.grp.values] + b[y_obs.grp.values] * y_obs.x
y = pm.Bernoulli('y', p=tt.nnet.sigmoid(est), observed = y_obs.y)
这是一个逻辑模型，而不是概率模型。如果你出于某种原因想要概率，你可以用标准概率函数替换 tt.nnet.sigmoid。
你的数据集仍然有点困难，但我认为这是因为数据生成中的一个错误：你假设所有 0.1 组的常数 a，但在模型中你允许 a 值因组而异。采样器在 sig_a 的值非常小时遇到问题（毕竟真正的值为 0...）。
编辑：更多解释：更改为使用标准法线 a_raw 和 b_raw 然后将它们转换为 Normal(mu=mu_a, sd=sig_a) 使用pm.Deterministic 不改变后验，但它使采样器更容易。它被称为非中心参数化。有关该主题的更深入描述，请参见例如 http://mc-stan.org/documentation/case-studies/divergences_and_bias.html，这也应该有助于您理解，为什么非常小的差异会产生问题。
编辑：生成新数据
n = 100
evts=np.random.randint(10,100,n)
x_g=np.random.binomial(1,.3,n)
i = np.zeros(evts.sum())
mu = np.random.normal(.5,.2,n)
mu0 = np.random.normal(.1,.05,n)
y_obs = pd.DataFrame({'y':i.copy(),'x':i.copy(),'grp':i.copy()},
index = range(evts.sum()))
i = 0
for grp in range(100):
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
if x_g[grp] ==1:
est = mu0[grp] + mu[grp]
else:
est=mu0[grp]
p_hat = tt.nnet.sigmoid(est).eval()
y_obs.loc[ind,'y_hat'] = est
y_obs.loc[ind,'y'] = np.random.binomial(1,p_hat,len(ind))
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs['grp']=y_obs.grp.astype(np.int64)
采样使用
with test_model:
trace = pm.sample(2000, tune=1000, njobs=4)
大约三分钟后完成
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
8%|▊ | 15977/200000 [00:15<02:51, 1070.66it/s]Median ELBO converged.
Finished [100%]: Average ELBO = -4,458.8
100%|██████████| 2000/2000 [02:48<00:00, 9.99it/s]
无分歧转换：
test_trace[1000:].diverging.sum()
全部使用pymc3和theano master。 （两者都准备好发布新版本）
我正在尝试通过添加概率转换并使结果服从伯努利分布，将层次模型从 Gelman 和 Hill reproduced in PyMC3 here 扩展到二进制结果数据。现在我正在使用玩具数据，所以我知道真正的价值。 Alpha 应该是 .1，beta 应该是 .5。
模型 运行 在扩展前使用 NUTS 采样器没问题。一旦我添加它，估计值就会缓慢增加并持续增加，直到模型在 10 到 200 次迭代之间停滞不前。这是它一直到 120（相对较长 运行）时的图像。
在扩展之前，Metropolis 需要进行 200,000 次迭代才能很好地修复真实参数值，但它最终做到了。扩展后它停滞在 30k 和 50k 之间。与 NUTS 不同，当您尝试停止它时它会完全崩溃 post-stall，所以我没有图片。更早地停止它给出了一个大致超过零的 beta 估计值，但分布范围很广。
代码贴在下面。
不知道是采样问题，还是规格问题。有没有更好的方法来指定 Probit？关于其他采样器的任何提示可以尝试吗？我已经尽可能地剥离我的模型以进行测试，并在我添加概率扩展后将其缩小到破坏，但我不知道下一步该怎么做。
#Generate Data
n=100
#Determine how many observations per group
evts=np.random.randint(10,100,n)
#Determine which groups will be receive treatment
x_g=np.random.binomial(1,.3,n)
#pre-create distribution of betas for groups
mu = np.random.normal(.5,.2,n)
#preallocate space in a dataframe
i = np.zeros(evts.sum())
y_obs = pd.DataFrame({'y':i.copy(),
'x':i.copy(),
'grp':i.copy()},
index = range(evts.sum()))
#populate dataframe with simulated data
i=0
for grp in range(100):
#index of observations for a given group
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
#generate outcomes using
#different dgp depending on treatment
if x_g[grp] ==1:
#shortcut to make sure 1>p>0
p_i = max((.1 + mu[grp]),0.01)
p_i = min(p_i,1)
out = np.random.binomial(1,p_i,evts[grp])
else:
out = np.random.binomial(1,.1,evts[grp])
#Assign to dataframe
y_obs.loc[ind,'y'] = out
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs = y_obs.astype(int)
print('starting model')
with pm.Model() as test_model:
#hyperpriors
mu_a=pm.Normal('mu_a',mu=0, sd=100**2)
sig_a = pm.Uniform('sig_a',lower=0,upper=100)
mu_b=pm.Normal('mu_b',mu=0, sd=100**2)
sig_b = pm.Uniform('sig_b',lower=0,upper=100)
#priors
a = pm.Normal('a',mu=mu_a,sd = sig_a, shape=n)
b = pm.Normal('b',mu=mu_b,sd = sig_b, shape=n)
eps = pm.Uniform('eps',lower=0,upper=100)
est = a[y_obs.grp] + b[y_obs.grp] * y_obs.x
#I get correct estimates when I
#stop here using commented out line.
# y_hat = pm.Normal('y_hat',
# mu=est,
# sd=eps,
# observed = y_obs.y)
#Probit transformation:
y_hat = pm.Normal('y_hat',
mu=est,
sd=eps,
shape=y_obs.shape[0])
mu_y = tt.mean(y_hat)
eps_hat = tt.var(y_hat)
p_hat = 0.5 * (1 + tt.erf((y_hat-mu_y) / (eps_hat*tt.sqrt(2))))
y = pm.Bernoulli('y',p=p_hat, observed = y_obs.y)
with test_model:
#Either:
mu,sds,elbo = pm.variational.advi(n=100000)
step = pm.NUTS(scaling=test_model.dict_to_array(sds),
is_cov=True)
test_trace = pm.sample(200, step, start=mu)
#or
# step=pm.Metropolis()
# test_trace = pm.sample(50000)
pm.traceplot(test_trace)#[-5000::3])
注意：编辑以修复行中的错字：'step = pm.NUTS(scaling=test_model.dict_to_array(sds),`
编辑：我为最初 posted 模型的概率扩展制作了更好的模拟数据。 （原始数据生成是现在 ADVI 给出了更好的估计，所以它从正确的位置开始，但 NUTS 仍然很快停止（大约十次迭代）。Metropolis 直接失败：我做了第一轮 5000 次迭代，并得到尝试绘制轨迹时出错。
新数据生成：
n=100
evts=np.random.randint(10,100,n)
x_g=np.random.binomial(1,.3,n)
i = np.zeros(evts.sum())
mu = np.random.normal(.5,.2,n)
mu0 = np.random.normal(.1,.05,n)
y_obs = pd.DataFrame({'y':i.copy(),'x':i.copy(),'grp':i.copy()},index = range(evts.sum()))
i=0
for grp in range(100):
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
if x_g[grp] ==1:
est = mu0[grp] + mu[grp]
else:
est=mu0[grp]
p_hat = tt.nnet.sigmoid(est).eval()
y_obs.loc[ind,'y_hat'] = est
y_obs.loc[ind,'y'] = np.random.binomial(1,p_hat,len(ind))
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs['grp']=y_obs.grp.astype(np.int64)
当 pymc3 尝试绘制密度时大都市出错：
ValueError: v cannot be empty
也许我误解了你的意思，但这个模型不应该工作：
with pm.Model() as test_model:
#hyperpriors
mu_a = pm.Flat('mu_a')
sig_a = pm.HalfCauchy('sig_a', beta=2.5)
mu_b = pm.Flat('mu_b')
sig_b = pm.HalfCauchy('sig_b', beta=2.5)
#priors
a_raw = pm.Normal('a_raw', mu=0, sd=1, shape=n)
a = pm.Deterministic('a', mu_a + sig_a * a_raw)
b_raw = pm.Normal('b_raw', mu=0, sd=1, shape=n)
b = pm.Deterministic('b', mu_b + sig_b * b_raw)
est = a[y_obs.grp.values] + b[y_obs.grp.values] * y_obs.x
y = pm.Bernoulli('y', p=tt.nnet.sigmoid(est), observed = y_obs.y)
这是一个逻辑模型，而不是概率模型。如果你出于某种原因想要概率，你可以用标准概率函数替换 tt.nnet.sigmoid。
你的数据集仍然有点困难，但我认为这是因为数据生成中的一个错误：你假设所有 0.1 组的常数 a，但在模型中你允许 a 值因组而异。采样器在 sig_a 的值非常小时遇到问题（毕竟真正的值为 0...）。
编辑：更多解释：更改为使用标准法线 a_raw 和 b_raw 然后将它们转换为 Normal(mu=mu_a, sd=sig_a) 使用pm.Deterministic 不改变后验，但它使采样器更容易。它被称为非中心参数化。有关该主题的更深入描述，请参见例如 http://mc-stan.org/documentation/case-studies/divergences_and_bias.html，这也应该有助于您理解，为什么非常小的差异会产生问题。
编辑：生成新数据
n = 100
evts=np.random.randint(10,100,n)
x_g=np.random.binomial(1,.3,n)
i = np.zeros(evts.sum())
mu = np.random.normal(.5,.2,n)
mu0 = np.random.normal(.1,.05,n)
y_obs = pd.DataFrame({'y':i.copy(),'x':i.copy(),'grp':i.copy()},
index = range(evts.sum()))
i = 0
for grp in range(100):
ind = list(range(i,(i+evts[grp])))
i += evts[grp]
if x_g[grp] ==1:
est = mu0[grp] + mu[grp]
else:
est=mu0[grp]
p_hat = tt.nnet.sigmoid(est).eval()
y_obs.loc[ind,'y_hat'] = est
y_obs.loc[ind,'y'] = np.random.binomial(1,p_hat,len(ind))
y_obs.loc[ind,'x'] = x_g[grp]
y_obs.loc[ind,'grp'] = grp
y_obs['grp']=y_obs.grp.astype(np.int64)
采样使用
with test_model:
trace = pm.sample(2000, tune=1000, njobs=4)
大约三分钟后完成
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
8%|▊ | 15977/200000 [00:15<02:51, 1070.66it/s]Median ELBO converged.
Finished [100%]: Average ELBO = -4,458.8
100%|██████████| 2000/2000 [02:48<00:00, 9.99it/s]
无分歧转换：
test_trace[1000:].diverging.sum()
全部使用pymc3和theano master。 （两者都准备好发布新版本）