pymc3 是否与 scikit-learn 中的 predict 方法等效？

如何解决pymc3 是否与 scikit-learn 中的 predict 方法等效？

我正在学习使用 PyMC3 的统计反思课程。在第 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("σ",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,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("α",a + b * shared_weights_s)
    sigma = pm.Uniform("σ",observed=adults.height)
    trace_adult = pm.sample()

现在，对于 shared_weights 的新值，我们有三个选择：

一件一件地做事情
替换为未知权重
附加到末尾

对于一对一的情况：

missing_weights = np.array([45,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,"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,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,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,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("α",a + b * weight_s)
    sigma = pm.Uniform("σ",observed=adults.height)
    trace_adult = pm.sample()

然后，在尝试时，我们pm.set_data()：

missing_weights = np.array([45,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

解决方法

我想我找到了答案：使用 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("σ",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

pymc3 是否与 scikit-learn 中的 predict 方法等效？

如何解决pymc3 是否与 scikit-learn 中的 predict 方法等效？

解决方法

相关推荐