pymc 对多个变量进行观察

Question

我正在使用黑客贝叶斯方法的线性回归示例，但无法将其扩展到我的用途。

我对一个随机变量进行了观察，对该随机变量进行了假设分布，最后对我观察到的该随机变量进行了另一个假设分布。a我如何尝试使用和上的中间分布对其进行建模b，但它抱怨Wrong number of dimensions: expected 0, got 1 with shape (788,).

为了描述实际模型，我预测一定数量 (n) 的培养电子邮件的转化率。我的先前观点是，转化率（由 和 上的 Beta 函数描述alpha）将通过某些因子 (0,inf]和beta进行更新alpha和缩放，这些因子从 n=0 的 1 开始，并在某个阈值处增加到最大值。betaab

# 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)

有办法让它发挥作用吗？

Answer 1

数据

首先，请注意，重复伯努利试验的总可能性正是二项式可能性，因此无需扩展到数据中的单个试验。我还建议使用 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]
})

解决方案

这将有助于简化模型，但解决方案实际上是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)

在这里，我们通过 sigmoid 将预测变量空间映射到概率空间，该 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)