1

我对 pymc3 完全陌生，所以请原谅这可能是微不足道的。我有一个非常简单的模型，我在其中预测二元响应函数。该模型几乎是此示例的逐字复制：https ://github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gelman_bioassay.py

我取回了模型参数（alpha、beta 和 theta），但我似乎无法弄清楚如何将模型的预测与输入数据重叠。我试过这样做（使用生物测定模型的说法）：

 from scipy.stats import binom

 mean_alpha = mean(trace['alpha'])
 mean_beta = mean(trace['beta'])

 pred_death = binom.rvs(n, 1./(1.+np.exp(-(mean_alpha + mean_beta * dose))))

然后绘制剂量与 pred_death 的关系，但这显然是不正确的，因为我每次都得到不同的二项分布图。

与此相关的是另一个问题，我如何评估拟合优度？在“入门” pymc3 教程中，我似乎找不到任何效果。

非常感谢您的任何建议！

4

2 回答 2

1

嗨，一个简单的方法如下：

from pymc3 import *
from numpy import ones, array

# Samples for each dose level
n = 5 * ones(4, dtype=int)
# Log-dose
dose = array([-.86, -.3, -.05, .73])


def invlogit(x):
    return np.exp(x) / (1 + np.exp(x))



with Model() as model:

    # Logit-linear model parameters
    alpha = Normal('alpha', 0, 0.01)
    beta = Normal('beta', 0, 0.01)

    # Calculate probabilities of death
    theta = Deterministic('theta', invlogit(alpha + beta * dose))

    # Data likelihood
    deaths = Binomial('deaths', n=n, p=theta, observed=[0, 1, 3, 5])
    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start, progressbar=True)




import matplotlib.pyplot as plt

death_fit = np.percentile(trace.theta,50,axis=0)

plt.plot(dose, death_fit,'g', marker='.', lw='1.25', ls='-', ms=5, mew=1)

plt.show()
于 2016-08-04T20:22:24.223 回答
0

如果要绘制剂量与 pred_death 的关系图，其中 pred_death 是根据 alpha 和 beta 的平均估计值计算得出的，则执行以下操作：

pred_death = 1./(1. + np.exp(-(mean_alpha + mean_beta * dose)))
plt.plot(dose, pred_death)

相反，如果您想绘制剂量与 pred_death 的关系，其中 pred_death 的计算考虑了 alpha 和 beta 的后验不确定性。那么可能最简单的方法是使用该功能sample_ppc：

可能是这样的

ppc = pm.sample_ppc（跟踪，样本=100，模型=pmmodel）

for i in range(100):
    plt.plot(dose, ppc['deaths'][i], 'bo', alpha=0.5)

使用后验预测检查 (ppc) 是一种通过将模型的预测与实际数据进行比较来检查模型行为的方法。这里有一个例子sample_ppc

其他选项可能是绘制平均值加上一些感兴趣的区间。

于 2016-05-21T00:34:50.237 回答