pymc 对多个变量的观察
pymc with observations on multiple variables
我正在使用 an example of linear regression from bayesian methods for hackers，但无法将其扩展到我的用途。
我有一个随机变量的观察值，一个关于该随机变量的假设分布，最后是我观察到的那个随机变量的另一个假设分布。我如何尝试使用 a 和 b 上的中间分布对其进行建模，但它抱怨 Wrong number of dimensions: expected 0, got 1 with shape (788,).
为了描述实际模型，我预测一定数量（n）的培养邮件的转化率。我的先验是转换率（由 alpha 和 beta 上的 Beta 函数描述）将通过 alpha 和 beta 按某些因素（0，inf）缩放来更新] a 和 b，当 n=0 时从 1 开始，并在某个阈值处增加到最大值。
# Generate predictive data, X and target data, Y
data = [
{'n': 0 , 'trials': 120, 'successes': 1},
{'n': 5 , 'trials': 111, 'successes': 2},
{'n': 10, 'trials': 78 , 'successes': 1},
{'n': 15, 'trials': 144, 'successes': 3},
{'n': 20, 'trials': 280, 'successes': 7},
{'n': 25, 'trials': 55 , 'successes': 1}]
X = np.empty(0)
Y = np.empty(0)
for dat in data:
X = np.insert(X, 0, np.ones(dat['trials']) * dat['n'])
target = np.zeros(dat['trials'])
target[:dat['successes']] = 1
Y = np.insert(Y, 0, target)
with pm.Model() as model:
alpha = pm.Uniform("alpha_n", 5, 13)
beta = pm.Uniform("beta_n", 1000, 1400)
n_sat = pm.Gamma("n_sat", alpha=20, beta=2, testval=10)
a_gamma = pm.Gamma("a_gamma", alpha=18, beta=15)
b_gamma = pm.Gamma("b_gamma", alpha=18, beta=27)
a_slope = pm.Deterministic('a_slope', 1 + (X/n_sat)*(a_gamma-1))
b_slope = pm.Deterministic('b_slope', 1 + (X/n_sat)*(b_gamma-1))
a = pm.math.switch(X >= n_sat, a_gamma, a_slope)
b = pm.math.switch(X >= n_sat, b_gamma, b_slope)
p = pm.Beta("p", alpha=alpha*a, beta=beta*b)
observed = pm.Bernoulli("observed", p, observed=Y)
有没有办法让它工作？
数据
首先，请注意重复伯努利试验的总似然是二项式似然，因此无需扩展到数据中的个别试验。我还建议使用 Pandas DataFrame 来管理您的数据 - 它有助于保持整洁：
import pandas as pd
df = pd.DataFrame({
'n': [0, 5, 10, 15, 20, 25],
'trials': [120, 111, 78, 144, 280, 55],
'successes': [1, 2, 1, 3, 7, 1]
})
解决方案
这将有助于简化模型，但真正的解决方案是向 p 随机变量添加一个 shape 参数，以便 PyMC3 知道如何解释一维参数。事实上，您确实希望每个 n 个案例都有不同的 p 分布，因此这里没有概念上的错误。
with pm.Model() as model:
# conversion rate hyperparameters
alpha = pm.Uniform("alpha_n", 5, 13)
beta = pm.Uniform("beta_n", 1000, 1400)
# switchpoint prior
n_sat = pm.Gamma("n_sat", alpha=20, beta=2, testval=10)
a_gamma = pm.Gamma("a_gamma", alpha=18, beta=15)
b_gamma = pm.Gamma("b_gamma", alpha=18, beta=27)
# NB: I removed pm.Deterministic b/c (a|b)_slope[0] is constant
# and this causes issues when using ArViZ
a_slope = 1 + (df.n.values/n_sat)*(a_gamma-1)
b_slope = 1 + (df.n.values/n_sat)*(b_gamma-1)
a = pm.math.switch(df.n.values >= n_sat, a_gamma, a_slope)
b = pm.math.switch(df.n.values >= n_sat, b_gamma, b_slope)
# conversion rates
p = pm.Beta("p", alpha=alpha*a, beta=beta*b, shape=len(df.n))
# observations
pm.Binomial("observed", n=df.trials, p=p, observed=df.successes)
trace = pm.sample(5000, tune=10000)
这个样本很好
并产生合理的转化率区间
但是 alpha_n 和 beta_n 的后验概率正好到达您之前的边界这一事实有点令人担忧：
我认为这样做的原因是，对于每个条件你只做 55-280 次试验，如果条件是独立的（最坏的情况），共轭会告诉我们你的 Beta 超参数应该在那个范围内.由于您正在进行回归，因此跨试验共享信息的最佳案例场景会将您的超参数置于试验总和 (788) 的范围内——但这是一个上限。因为你超出了这个范围，这里的问题是你强迫模型的估计比你真正有证据支持的更精确。但是，如果先验是基于强有力的 independent 证据，则可以证明这一点。
否则，我建议扩大影响最终 alpha*a 和 beta*b 数字的那些先验的范围（这些数字的总和应该接近你在后验中的试验计数）。
替代模型
我可能会按照以下几行做一些事情，我认为它具有更透明的参数化，尽管它与您的模型并不完全相同：
with pm.Model() as model_br_sp:
# regression coefficients
alpha = pm.Normal("alpha", mu=0, sd=1)
beta = pm.Normal("beta", mu=0, sd=1)
# saturation parameters
saturation_point = pm.Gamma("saturation_point", alpha=20, beta=2)
max_success_rate = pm.Beta("max_success_rate", 1, 9)
# probability of conversion
success_rate = pm.Deterministic("success_rate",
pm.math.switch(df.n.values > saturation_point,
max_success_rate,
max_success_rate*pm.math.sigmoid(alpha + beta*df.n)))
# observations
pm.Binomial("successes", n=df.trials, p=success_rate, observed=df.successes)
trace_br_sp = pm.sample(draws=5000, tune=10000)
在这里，我们通过以最大成功率最大化的 sigmoid 将预测变量 space 映射到概率 space。饱和点的先验与您的相同，而最大成功率的先验提供的信息很少（Beta [1,9] - 尽管我会说它几乎也以平坦的先验运行）。这也很好采样，
并给出了相似的间隔（尽管切换点似乎占主导地位）：
我们可以比较两个模型，发现它们的解释力没有显着差异：
import arviz as az
model_compare = az.compare({'Binomial Regression w/ Switchpoint': trace_br_sp,
'Original Model': trace})
az.plot_compare(model_compare)
我正在使用 an example of linear regression from bayesian methods for hackers，但无法将其扩展到我的用途。
我有一个随机变量的观察值，一个关于该随机变量的假设分布，最后是我观察到的那个随机变量的另一个假设分布。我如何尝试使用 a 和 b 上的中间分布对其进行建模，但它抱怨 Wrong number of dimensions: expected 0, got 1 with shape (788,).
为了描述实际模型，我预测一定数量（n）的培养邮件的转化率。我的先验是转换率（由 alpha 和 beta 上的 Beta 函数描述）将通过 alpha 和 beta 按某些因素（0，inf）缩放来更新] a 和 b，当 n=0 时从 1 开始，并在某个阈值处增加到最大值。
# Generate predictive data, X and target data, Y
data = [
{'n': 0 , 'trials': 120, 'successes': 1},
{'n': 5 , 'trials': 111, 'successes': 2},
{'n': 10, 'trials': 78 , 'successes': 1},
{'n': 15, 'trials': 144, 'successes': 3},
{'n': 20, 'trials': 280, 'successes': 7},
{'n': 25, 'trials': 55 , 'successes': 1}]
X = np.empty(0)
Y = np.empty(0)
for dat in data:
X = np.insert(X, 0, np.ones(dat['trials']) * dat['n'])
target = np.zeros(dat['trials'])
target[:dat['successes']] = 1
Y = np.insert(Y, 0, target)
with pm.Model() as model:
alpha = pm.Uniform("alpha_n", 5, 13)
beta = pm.Uniform("beta_n", 1000, 1400)
n_sat = pm.Gamma("n_sat", alpha=20, beta=2, testval=10)
a_gamma = pm.Gamma("a_gamma", alpha=18, beta=15)
b_gamma = pm.Gamma("b_gamma", alpha=18, beta=27)
a_slope = pm.Deterministic('a_slope', 1 + (X/n_sat)*(a_gamma-1))
b_slope = pm.Deterministic('b_slope', 1 + (X/n_sat)*(b_gamma-1))
a = pm.math.switch(X >= n_sat, a_gamma, a_slope)
b = pm.math.switch(X >= n_sat, b_gamma, b_slope)
p = pm.Beta("p", alpha=alpha*a, beta=beta*b)
observed = pm.Bernoulli("observed", p, observed=Y)
有没有办法让它工作？
数据
首先，请注意重复伯努利试验的总似然是二项式似然，因此无需扩展到数据中的个别试验。我还建议使用 Pandas DataFrame 来管理您的数据 - 它有助于保持整洁：
import pandas as pd
df = pd.DataFrame({
'n': [0, 5, 10, 15, 20, 25],
'trials': [120, 111, 78, 144, 280, 55],
'successes': [1, 2, 1, 3, 7, 1]
})
解决方案
这将有助于简化模型，但真正的解决方案是向 p 随机变量添加一个 shape 参数，以便 PyMC3 知道如何解释一维参数。事实上，您确实希望每个 n 个案例都有不同的 p 分布，因此这里没有概念上的错误。
with pm.Model() as model:
# conversion rate hyperparameters
alpha = pm.Uniform("alpha_n", 5, 13)
beta = pm.Uniform("beta_n", 1000, 1400)
# switchpoint prior
n_sat = pm.Gamma("n_sat", alpha=20, beta=2, testval=10)
a_gamma = pm.Gamma("a_gamma", alpha=18, beta=15)
b_gamma = pm.Gamma("b_gamma", alpha=18, beta=27)
# NB: I removed pm.Deterministic b/c (a|b)_slope[0] is constant
# and this causes issues when using ArViZ
a_slope = 1 + (df.n.values/n_sat)*(a_gamma-1)
b_slope = 1 + (df.n.values/n_sat)*(b_gamma-1)
a = pm.math.switch(df.n.values >= n_sat, a_gamma, a_slope)
b = pm.math.switch(df.n.values >= n_sat, b_gamma, b_slope)
# conversion rates
p = pm.Beta("p", alpha=alpha*a, beta=beta*b, shape=len(df.n))
# observations
pm.Binomial("observed", n=df.trials, p=p, observed=df.successes)
trace = pm.sample(5000, tune=10000)
这个样本很好
并产生合理的转化率区间
但是 alpha_n 和 beta_n 的后验概率正好到达您之前的边界这一事实有点令人担忧：
我认为这样做的原因是，对于每个条件你只做 55-280 次试验，如果条件是独立的（最坏的情况），共轭会告诉我们你的 Beta 超参数应该在那个范围内.由于您正在进行回归，因此跨试验共享信息的最佳案例场景会将您的超参数置于试验总和 (788) 的范围内——但这是一个上限。因为你超出了这个范围，这里的问题是你强迫模型的估计比你真正有证据支持的更精确。但是，如果先验是基于强有力的 independent 证据，则可以证明这一点。
否则，我建议扩大影响最终 alpha*a 和 beta*b 数字的那些先验的范围（这些数字的总和应该接近你在后验中的试验计数）。
替代模型
我可能会按照以下几行做一些事情，我认为它具有更透明的参数化，尽管它与您的模型并不完全相同：
with pm.Model() as model_br_sp:
# regression coefficients
alpha = pm.Normal("alpha", mu=0, sd=1)
beta = pm.Normal("beta", mu=0, sd=1)
# saturation parameters
saturation_point = pm.Gamma("saturation_point", alpha=20, beta=2)
max_success_rate = pm.Beta("max_success_rate", 1, 9)
# probability of conversion
success_rate = pm.Deterministic("success_rate",
pm.math.switch(df.n.values > saturation_point,
max_success_rate,
max_success_rate*pm.math.sigmoid(alpha + beta*df.n)))
# observations
pm.Binomial("successes", n=df.trials, p=success_rate, observed=df.successes)
trace_br_sp = pm.sample(draws=5000, tune=10000)
在这里，我们通过以最大成功率最大化的 sigmoid 将预测变量 space 映射到概率 space。饱和点的先验与您的相同，而最大成功率的先验提供的信息很少（Beta [1,9] - 尽管我会说它几乎也以平坦的先验运行）。这也很好采样，
并给出了相似的间隔（尽管切换点似乎占主导地位）：
我们可以比较两个模型，发现它们的解释力没有显着差异：
import arviz as az
model_compare = az.compare({'Binomial Regression w/ Switchpoint': trace_br_sp,
'Original Model': trace})
az.plot_compare(model_compare)