我使用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一个测试功能所描述的）将通过具有alpha和beta由一些因素（0，缩放被更新INF] a和b，其中从1开始对n = 0，并增加至在某个阈值的最大值。
# 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)
有没有办法得到这个工作？
首先，注意，重复的伯努利试验的总可能性完全是一个二项式可能性，所以没有必要扩大到单个试验中的数据。我也建议使用熊猫数据帧来管理你的数据 - 这是有助于保持整洁：
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]
})
这将有助于简化模型，但解决的办法确实是一个shape参数添加到p随机变量，以便PyMC3知道如何解释的一组维的参数。事实是，你要为你的每个p有不同的n分布，所以没有什么概念错在这里。
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)
在这里，我们通过在最大成功率马克塞斯S形映射预测空间概率空间。在饱和点的现有等同于你的，而对最大的成功率是弱信息（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)