我正在使用an example of linear regression from bayesian methods for hackers但无法将其扩展到我的用途。
我对一个随机变量进行了观察,对该随机变量进行了假设分布,最后对我观察到的该随机变量进行了另一个假设分布。我如何尝试使用 a
和 b
上的中间分布对其进行建模,但它提示维数错误:预期为 0,形状为 1 (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) 的范围内 - 但这是上限。因为您超出了这个范围,所以这里的问题是您迫使模型的估计比您真正有证据支持的更精确。然而,如果先验基于强有力的独立证据,就可以证明这一点是合理的。
否则,我建议扩大影响最终 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 将预测空间映射到概率空间。饱和点的先验与你的相同,而最大成功率的先验信息很少(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)
关于python - pymc 对多个变量进行观察,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/54545661/