我正在学习使用 PyMC3 的 Statistical Rethinking 课程。在第 4 章结束时，他们要求原始 (!Kung) 数据集中没有的数据点的单个值的 HDI。可以在 PyMC3 中做到这一点吗？
在 scikit-learn 中，您拥有fit()和predict()，并且您能够预测全新输入的输出。
使用 PyMC3，您sample()可以获得跟踪，并且可以要求进行后验预测检查，但我无法将任何参数传递给我感兴趣的值。我确实设法以迂回的方式做到了使用共享 theano 变量，也可以手动使用。
编辑：我在最后添加了一个pm.Data()示例pm.set_data()。我认为这可能是答案，但在将其标记为已回答之前，我正在等待其他人确认。
这就是我所做的。
weight_s是标准化的重量数据。标准化是通过这个函数完成的：
def standardize(array, reference=None):
if reference is None:
reference = array
return (array - reference.mean()) / reference.std()
这是 PyMC3 模型：
import pymc3 as pm
with pm.Model() as m_adult:
a = pm.Normal("α", mu=155, sd=20)
b = pm.Lognormal("β", mu=0, sd=1)
mu = pm.Deterministic("μ", a + b * adults.weight_s)
sigma = pm.Uniform("σ", 0, 50)
height = pm.Normal("height", mu=mu, sd=sigma, observed=adults.height)
trace_adult = pm.sample()
这是数据的样子（注意 HDI 是 89%）：
height_pred = pm.fast_sample_posterior_predictive(trace_adult, model=m_adult)["height"]
fig, ax = plt.subplots()
ax.plot(adults.weight, adults.height, ".")
ax.plot(adults.weight, trace_adult.μ.mean(axis=0), color="black")
az.plot_hdi(adults.weight, trace_adult.μ, ax=ax, color="black")
az.plot_hdi(adults.weight, height_pred, ax=ax)
ax.set(xlabel="weight", ylabel="height")
fig.tight_layout()
首先，我将向您展示手动版本：
missing_weights = np.array([45, 40, 65, 31, 53])
expected_height = np.array([
(trace_adult.α + trace_adult.β * standardize(weight, adults.weight)).mean()
for weight in missing_weights
])
hdis = np.array([
az.hdi(np.random.normal(
trace_adult.α + trace_adult.β * standardize(weight, adults.weight),
trace_adult.σ,
)) for weight in missing_weights
])
data = np.vstack((missing_weights, expected_height, hdis.T)).T
missing_df = pd.DataFrame(, columns=["weight", "expected_height", "hdi_lower", "hdi_upper"])
print(missing_df)
这给了我们：
weight expected_height hdi_lower hdi_upper
0 45 154.603176 146.981285 163.149938
1 40 150.105295 142.095583 158.474277
2 65 172.594698 164.401102 180.786641
3 31 142.009110 134.163952 150.233028
4 53 161.799785 153.881956 170.209779
如果您查看图表，这些数字是有意义的。
现在为共享变量。我们可以这样修改模型：
from theano import shared
shared_weights_s = shared(adults.weight_s.values)
with pm.Model() as m_adult:
a = pm.Normal("α", mu=155, sd=20)
b = pm.Lognormal("β", mu=0, sd=1)
mu = pm.Deterministic("μ", a + b * shared_weights_s)
sigma = pm.Uniform("σ", 0, 50)
height = pm.Normal("height", mu=mu, sd=sigma, observed=adults.height)
trace_adult = pm.sample()
现在，对于 shared_weights 的新值，我们有三个选择：
- 一件一件地做事
- 替换为未知权重
- 追加到最后
对于一对一的情况：
missing_weights = np.array([45, 40, 65, 31, 53])
rows = []
for weight in missing_weights:
row = [weight]
shared_weights_s.set_value(standardize(np.array([weight]), adults.weight))
height_pred_single = pm.fast_sample_posterior_predictive(trace_adult, model=m_adult)["height"]
row.append(height_pred_single.mean())
row.extend(list(az.hdi(height_pred_single).mean(axis=0)))
rows.append(row)
missing_df = pd.DataFrame(rows, columns=["weight", "expected_height", "hdi_lower", "hdi_upper"])
print(missing_df)
为他们所有人做这件事给了我们：
weight expected_height hdi_lower hdi_upper
0 45 154.604520 146.485327 162.713345
1 40 150.113378 142.001151 158.263953
2 65 172.580212 164.357970 180.843184
3 31 142.010954 133.786200 150.142080
4 53 161.792962 153.651266 169.926615
您可以一次完成所有操作：
missing_weights = np.array([45, 40, 65, 31, 53])
shared_weights_s.set_value(standardize(missing_weights, adults.weight))
height_pred_replace = pm.fast_sample_posterior_predictive(trace_adult, model=m_adult)["height"]
missing_df = pd.DataFrame(missing_weights, columns=["weight"])
missing_df["expected_height"] = height_pred_replace.mean(axis=0)
missing_df[["hdi_lower", "hdi_upper"]] = az.hdi(height_pred_replace)
print(missing_df)
这给了我们：
weight expected_height hdi_lower hdi_upper
0 45 154.578096 147.066342 163.069805
1 40 150.042506 141.561599 158.120596
2 65 172.568430 164.079591 180.536870
3 31 142.080048 134.173959 150.345556
4 53 161.830472 153.327694 169.717058
最后，我们可以将它添加到之前共享权重变量的末尾并取尾部：
missing_weights = np.array([45, 40, 65, 31, 53])
shared_weights_s.set_value(np.append(adults.weight_s.values, standardize(missing_weights, adults.weight)))
height_pred_append = pm.fast_sample_posterior_predictive(trace_adult, model=m_adult)["height"]
missing_df = pd.DataFrame(missing_weights, columns=["weight"])
missing_df["expected_height"] = height_pred_append.mean(axis=0)[-len(missing_weights):]
missing_df[["hdi_lower", "hdi_upper"]] = az.hdi(height_pred_append)[-len(missing_weights):]
print(missing_df)
这给了我们：
weight expected_height hdi_lower hdi_upper
0 45 154.640287 146.093825 162.477313
1 40 150.088713 142.168331 158.314038
2 65 172.633776 164.086280 180.483805
3 31 142.019331 133.516545 150.491937
4 53 161.880175 153.530868 169.771088
如您所见，所有这些方法最终都会给出相同的结果。有没有官方/最好的方法？可以不设置全局共享变量并修改它吗？PyMC3有这样的功能吗，还是以后可能会添加的东西？（如果足够简单，我也许可以为此提出拉取请求；我还是 PyMC3 的新手。）
编辑：我想我找到了答案：使用pm.Data().
with pm.Model() as m_adult:
weight_s = pm.Data("weight_s", adults.weight_s.values)
a = pm.Normal("α", mu=155, sd=20)
b = pm.Lognormal("β", mu=0, sd=1)
mu = pm.Deterministic("μ", a + b * weight_s)
sigma = pm.Uniform("σ", 0, 50)
height = pm.Normal("height", mu=mu, sd=sigma, observed=adults.height)
trace_adult = pm.sample()
然后，在尝试时，我们pm.set_data()：
missing_weights = np.array([45, 40, 65, 31, 53])
with m_adult:
pm.set_data({"weight_s": standardize(missing_weights, adults.weight)})
height_pred_data = pm.fast_sample_posterior_predictive(trace_adult)["height"]
missing_df = pd.DataFrame(missing_weights, columns=["weight"])
missing_df["expected_height"] = height_pred_data.mean(axis=0)
missing_df[["hdi_lower", "hdi_upper"]] = az.hdi(height_pred_data)
print(missing_df)
这使：
weight expected_height hdi_lower hdi_upper
0 45 154.584063 145.828088 162.512174
1 40 150.184853 142.272258 158.451555
2 65 172.662069 164.522903 180.803430
3 31 141.949137 133.310865 149.811098
4 53 161.719867 153.848599 169.638495