TowardsDataScience 博客中文翻译 2021（一百三十四）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

贝叶斯 AB 检验—第二部分—收入

原文：https://towardsdatascience.com/bayesian-ab-testing-part-ii-revenue-1fbcf04f96cd?source=collection_archive---------10-----------------------

米歇尔·汉德森在 Unsplash 上的照片

系列结构

这篇文章是关于将贝叶斯 AB 测试方法应用于现实生活产品场景的系列博客文章的第二部分。它使用了在系列的第一部分中讨论的一些概念。

1. 基于转换的测试指标(速率指标)的建模和分析
2. 基于收入的测试指标的建模和分析(连续指标)
3. 计算测试持续时间
4. 选择一个合适的先验
5. 用多个变量运行测试

让我们开始讨论如何在贝叶斯产品实验中建模、使用和分析收入测试指标。

实验背景

假设我们最近对应用程序中的商店功能进行了 UX 更改。我们相信这些变化让我们的用户更容易进行更大的应用内购买，我们希望在向更广泛的用户群发布之前进行 AB 测试。我们假设我们所做的改变将导致每用户平均收入的显著提高。

我们将每个用户产生的收入建模为随机变量𝑅=𝑋∗𝑌，其中:

𝑋是一个伯努利随机变量，它指的是用户是否进行了购买，转换概率为𝜆:𝑋𝐵𝑒𝑟(𝜆)

𝑌是一个指数随机变量，指的是购买的规模，如果它是与利率参数𝜃:𝑌𝐸𝑥𝑝(𝜃)

在这个模型中，我们知道每笔销售的平均收入是由 1/𝜃给出的，每用户的平均收入是由𝜆/𝜃.给出的

我们分析的第一步是查看过去的数据，为该模型的关键参数𝜆和𝜃选择合适的先验分布。

接下来，我们设置一个损失阈值——在我们错误地推出错误变体的情况下，我们愿意接受的最大预期损失。与任何统计建模一样，贝叶斯实验方法是建立在真实世界数据的近似值之上的。因此，我们总是有可能从测试中得出错误的结论。这个损失阈值允许我们说，即使我们得出了错误的结论，我们也可以确信每个用户的平均收入不会下降超过这个阈值。

最后，我们以随机实验的形式抽取样本，并使用这些样本来更新分布，以及我们对商店控制和治疗版本下的𝜆和𝜃的信念。然后我们可以用这些后验分布来计算我们感兴趣的概率和预期损失。

为了给我们的参数选择先验分布，让我们先看看我们最近收集的一些过去的购买数据。我生成了一个样本先验数据集，我们可以用它来做这个练习。

import pandas as pd
import numpy as npprior_data = pd.read_csv('prior_data.csv')print(prior_data.head())
print(prior_data.shape)

因为这个数据集是人工生成的，所以它已经是这个练习的理想格式。在现实世界中，我们可能需要执行一些 ETL 操作来获得这种格式的数据。然而，这超出了本文的范围。

我们看到，我们有一个 5268 个用户的样本大小，对于每个用户，我们可以看到他们是否进行了购买。对于那些确实进行了购买的用户，我们还可以看到购买的规模。我们可以继续计算之前的转换率、每笔销售的平均收入、每笔销售收入的比率参数以及每用户的平均收入。

conversion_rate = prior_data['converted'].sum()/prior_data.shape[0]converted = prior_data[prior_data['converted'] == 1]
avg_purchase = converted['revenue'].mean()print(f'Prior Conversion Rate is {round(conversion_rate, 3)}. Average Revenue per Sale is {round(avg_purchase, 3)}.')print(f'Rate Parameter for Revenue per Sale is {round(1/avg_purchase, 3)}. Average Revenue per User is {round(conversion_rate*avg_purchase, 3)}.')

选择先验分布

利用上面的信息，我们可以选择𝜆和𝜃.的先验分布使用与前一篇文章相同的逻辑，我们可以为我们的转换概率𝜆.选择𝐵𝑒𝑡𝑎(7,15 的先验总之，我们选择使用贝塔分布，因为它是[0，1]上的灵活分布，并且在之前是一个好的共轭。

对于𝜃，每笔销售收入的比率参数，我们将使用伽马分布，因为它是一个在[0，∞)上的灵活分布，也是一个很好的共轭先验。当我们用实验数据计算后验概率时，这将使我们的计算更容易。

我们可以选择非常弱的先验分布𝐺𝑎𝑚𝑚𝑎(0.1，0.1)。[1]

设置损失阈值

既然我们已经选择了先验，我们需要选择我们的𝜖，这是我们在错误选择错误变量的情况下愿意接受的最高预期损失。让我们假设商店不是我们的主要收入来源，但它对我们非常重要，所以我们希望在这个𝜖.上保守一些我们可以设定𝜖 = 0.005。

我们有先验分布和预期损失的阈值，所以我们可以开始运行我们的实验并从中收集数据。

实验结果

让我们假设我们已经让我们的实验运行了几个星期，并想检查我们是否能从中得出任何结论。为了做到这一点，我们需要使用我们的实验数据来计算我们的后验分布，然后我们可以使用它来计算每个变量变得更好的概率，以及错误选择每个变量的预期损失。

出于本练习的目的，我生成了一个样本实验数据集。让我们从探索它开始。

experiment_data = pd.read_csv('experiment_data.csv')print(experiment_data.head())
print(experiment_data.shape)

我们看到该数据集与之前的数据集相似，只是多了一列，用于指定用户被分配到哪个组，以及他们看到哪个变量。再次值得注意的是，由于这个数据集是人工生成的，它已经是这个练习的理想格式，不需要额外的 ETL 操作。

我们现在可以开始汇总数据了。

results = experiment_data.groupby('group').agg({'userId': pd.Series.nunique, 'converted': sum, 'revenue': sum})results.rename({'userId': 'sampleSize'}, axis=1, inplace=True)results['conversionRate'] = results['converted']/results['sampleSize']results['revenuePerSale'] = results['revenue']/results['converted']
print(results)

经检查，似乎两组有相似的转换率，但治疗有更好的每销售收入。然而，为了更新我们关于𝜆和𝜃的两种变体的信念，我们需要进行进一步的计算。

使用前一篇文章的中概述的计算，我们可以计算两个变量的𝜆后验分布。

from scipy.stats import beta, gamma
import seaborn as sns
import matplotlib.pyplot as pltcontrol_cr = beta(7 + results.loc['control', 'converted'], 15 + results.loc['control', 'sampleSize'] - results.loc['control', 'converted'])treatment_cr = beta(7 + results.loc['treatment', 'converted'], 15 + results.loc['treatment', 'sampleSize'] - results.loc['treatment', 'converted'])fig, ax = plt.subplots()x = np.linspace(0,1,1000)ax.plot(x, control_cr.pdf(x), label='control')
ax.plot(x, treatment_cr.pdf(x), label='treatment')
ax.set_xlabel('Conversion Probability')
ax.set_ylabel('Density')
ax.set_title('Experiment Posteriors')
ax.legend()

转换后验概率(图片由作者提供)

𝜆_𝑐和𝜆_𝑡的后验分布几乎相同。通过观察我们可以看出，处理对转换概率没有大的影响。

让我们来看看它是如何影响收入的比率参数的。我们将使用下面的结果[2]来计算𝜃_𝑐和𝜃_𝑡:

假设先验

𝜃 ∼ 𝐺𝑎𝑚𝑚𝑎(𝑘, Θ)

假设向𝑛游客展示了一个变种。𝑐换算成每销售𝑠.的平均收入然后，后验分布由下式给出

𝜃|𝑐,𝑠 ∼ 𝐺𝑎𝑚𝑚𝑎(𝑘 + 𝑐, Θ/(1 + Θ𝑐𝑠))

(如果你对这个结果背后的数学原理感到好奇，这里有的证明。)

让我们继续计算我们的后验概率𝜃_𝑐和𝜃_𝑡.

control_rr = gamma(a=(0.1 + results.loc['control', 'converted']), scale=(0.1/(1 + (0.1)*results.loc['control', 'converted']*results.loc['control', 'revenuePerSale'])))treatment_rr = gamma(a=(0.1 + results.loc['treatment', 'converted']), scale=(0.1/(1 + (0.1)*results.loc['treatment', 'converted']*results.loc['treatment', 'revenuePerSale'])))fig, ax = plt.subplots()x = np.linspace(0,3,1000)ax.plot(x, control_rr.pdf(x), label='control')
ax.plot(x, treatment_rr.pdf(x), label='treatment')
ax.set_xlabel('Rate Parameter')
ax.set_ylabel('Density')
ax.set_title('Experiment Posteriors')
ax.legend()

费率参数(每笔销售收入)后验概率(图片由作者提供)

后验分布根本没有重叠，所以我们可以相当确定其中一个变量比另一个更好。处理率参数低于控制率参数，因此更好，因为它将导致更高的平均销售收入。为了更清楚地看到这一点，让我们画出在这些后验概率下每笔销售的平均收入的分布。

fig, ax = plt.subplots()x = np.linspace(0,3,1000)
z = [1/i for i in x]ax.plot(x, control_rr.pdf(z), label='control')
ax.plot(x, treatment_rr.pdf(z), label='treatment')
ax.set_xlabel('Avg Revenue per Sale')
ax.set_ylabel('Density')
ax.set_title('Experiment Posteriors')
ax.legend()

每笔销售的平均收入(图片由作者提供)

我们可以清楚地看到，治疗组的平均销售收入高于对照组。

所以两种变体有相似的转换率，但是 treatment 有更好的每笔销售的平均收入，所以它也必须有更好的每用户平均收入。让我们画出每个用户平均收入的后验概率来证实这一点。我们可以进行一些模拟来帮助我们。

control_conversion_simulation = np.random.beta(7 + results.loc['control', 'converted'], 15 + results.loc['control', 'sampleSize'] - results.loc['control', 'converted'], size=100000)treatment_conversion_simulation = np.random.beta(7 + results.loc['treatment', 'converted'], 15 + results.loc['treatment', 'sampleSize'] - results.loc['treatment', 'converted'], size=100000)control_revenue_simulation = np.random.gamma(shape=(0.1 + results.loc['control', 'converted']), scale=(0.1/(1 + (0.1)*results.loc['control', 'converted']*results.loc['control', 'revenuePerSale'])), size=100000)treatment_revenue_simulation = np.random.gamma(shape=(0.1 + results.loc['treatment', 'converted']), scale=(0.1/(1 + (0.1)*results.loc['treatment', 'converted']*results.loc['treatment', 'revenuePerSale'])), size=100000)control_avg_purchase = [i/j for i,j in zip(control_conversion_simulation, control_revenue_simulation)]treatment_avg_purchase = [i/j for i,j in zip(treatment_conversion_simulation, treatment_revenue_simulation)]fig, ax = plt.subplots()x = np.linspace(0,1,1000)ax.hist(control_avg_purchase, density=True, label='control', histtype='stepfilled', bins=100)
ax.hist(treatment_avg_purchase, density=True, label='treatment', histtype='stepfilled', bins=100)
ax.set_xlabel('Avg Revenue per User')
ax.set_ylabel('Density')
ax.set_title('Experiment Posteriors')
ax.legend()

每用户平均收入(图片由作者提供)

治疗组的平均收入明显高于对照组。通过检查这些后验，我们可以很有把握地认为治疗是更好的。然而，为了量化这一点，我们需要计算𝑝(𝜆_𝑡/𝜃_𝑡 ≥ 𝜆_𝑐/𝜃_𝑐)和𝐸𝐿，错误选择治疗的预期损失。

treatment_won = [i <= j for i,j in zip(control_avg_purchase, treatment_avg_purchase)]chance_to_beat_ctrl = np.mean(treatment_won)
print(f'Chance of beating control: {round(chance_to_beat_ctrl, 3)}.')

从模拟中我们看到𝑝(𝜆_𝑡/𝜃_𝑡 ≥ 𝜆_𝑐/𝜃_𝑐) = 1，所以治疗有 100%的机会比对照组好。

既然我们已经计算了治疗效果更好的可能性，我们需要计算𝐸𝐿.每个变量的损失函数由下式给出

我们用这个来计算预期损失。

loss_control = [max(j - i, 0) for i,j in zip(control_avg_purchase, treatment_avg_purchase)]loss_treatment = [max(i - j, 0) for i,j in zip(control_avg_purchase, treatment_avg_purchase)]all_loss_control = [int(i)*j for i,j in zip(treatment_won, loss_control)]all_loss_treatment = [(1 - int(i))*j for i,j in zip(treatment_won, loss_treatment)]expected_loss_control = np.mean(all_loss_control)
expected_loss_treatment = np.mean(all_loss_treatment)print(f'Expected loss of choosing control: {round(expected_loss_control, 3)}. Expected loss of choosing treatment: {round(expected_loss_treatment, 3)}')

从运行模拟中我们看到

𝐸𝐿=0 < 0.005=𝜖.

因为其中一个变异的预期损失低于我们在测试开始时设定的阈值，所以测试达到了显著性。我们可以很有把握地得出结论，治疗效果更好，错误选择治疗的预期成本不会超过我们能接受的水平。因此，我们强烈建议将追加销售屏幕的处理方式推广到我们的其他用户群。

我希望这个案例研究有助于您理解实现贝叶斯 AB 测试方法所需的计算。观看这个空间的系列的下一部分！

参考

[1] 通过 Omri Fima 的贝叶斯 A/B 测试优化收入

[2] VWO 白皮书作者:C.Stucchio

我这篇文章的代码可以在这里找到

感谢阅读这篇文章！我希望它能帮助您更好地理解如何为收入和其他计数指标实现贝叶斯 AB 测试方法。

如果你喜欢阅读我的文章，愿意支持我的写作，并且正在考虑订阅一个媒体，请随时使用我下面的推荐链接。我会从你的订阅费中提成。

https://medium.com/@kaushsk12/membership

贝叶斯 AB 测试第三部分测试持续时间

原文：https://towardsdatascience.com/bayesian-ab-testing-part-iii-test-duration-f2305215009c?source=collection_archive---------22-----------------------

在 Unsplash 上由 Aron 视觉拍摄的照片

系列结构

这篇文章是将贝叶斯 AB 测试方法应用于现实生活产品场景的系列博客文章的第三部分。它使用了在系列的第一部分中讨论的一些概念。

1. 基于转换的测试指标(速率指标)的建模和分析
2. 基于收入的测试指标(持续指标)的建模和分析
3. 计算测试持续时间
4. 选择一个合适的先验
5. 用多个变量运行测试

实验背景

根据上一篇文章中使用的示例，我们假设我们最近更改了追加销售屏幕上的信息，并希望在向更广泛的用户群发布之前进行 AB 测试。我们假设我们所做的改变将会带来更好的转化率。

我们继续将转换率建模为具有转换概率𝜆的伯努利随机变量，进而我们使用𝐵𝑒𝑡𝑎(7,15).的先验分布进行建模然后我们选择我们的预期损失阈值𝜖 = 0.0015。我们现在准备运行我们的测试，但是我们要运行多长时间呢？

在我们深入分析来回答这个问题之前，让我们首先考虑一下测试持续时间实际上意味着什么。以下公式给出了计算测试持续时间的常用方法:

测试中包含的变量数量和用户比例是在测试设计过程中就已经决定的因素。我们预计一周内的用户数量是基于前几周的数据进行的简单计算。所以本质上，问题归结为在实验中为每个变体选择一个合理的所需样本量。

人们经常建议，贝叶斯产品实验应该运行，直到其中一个变量的预期损失低于我们的阈值𝜖，在这一点上，我们宣布该变量为赢家。这样我们就不需要担心计算所需的样本量。然而，遵循这种方法可能会导致我们错误地选择一个变体来推出，因为这是一个被称为窥视 [1]的概念。为了进一步探讨这一点，让我们考虑下面的例子。

偷看

让我们考虑一下𝜆_𝑐和𝜆_𝑡非常相似的情况。我们模拟运行 AB 测试，并期望结果是非决定性的。然而，我们看到以下结果。

偷看例子(图片由作者提供)

我们从右侧的转换率中可以看出，𝜆_𝑐和𝜆_𝑡在前 1000 个样本中趋于一致。然而，更有趣的观察是关于左边的图表。我们看到，在前几个样本中，控制预期损失低于阈值。如果我们在这里停止测试，我们会得出控制明显优于治疗的结论，这将是错误的结论。我们还看到，如果我们在这里不停止测试，根据我们停止测试的时间，所有三种结果(对照组获胜、治疗组获胜和非决定性测试)都是可能的。那么我们如何决定何时停止测试呢？

运行我们实验的模拟将帮助我们避免落入窥视的陷阱，并最终帮助我们为每个变量选择一个合理的所需样本量。在我们研究这是如何工作的之前，让我们考虑下面的警告。

警告

下面提出的方法的主要警告是，它没有考虑任何选择测试持续时间的实验设计参数。特别是，它没有考虑任何季节或基于时间的转换变化。在现实世界中，您用作测试指标的转换很可能会根据一周中的某一天甚至一天中的某个时间而变化。因此，为了避免由于季节变化而得出错误的结论，至少运行几周的测试是值得的。我建议使用这篇文章中概述的方法来计算运行测试需要多长时间。

样本量计算

影响所需样本量的主要因素有三个:

• 可察觉的最小影响𝛿—𝜆的最小变化我们需要推广治疗
• 预期损失阈值𝜖——在我们错误宣布赢家的情况下，我们愿意接受的转化率的最大损失
• 转换概率𝜆标度

我们将在后面的文章中深入探讨这些因素是如何影响样本量的，但是让我们先建立一个计算合理样本量的过程。

让我们先给上述因素赋值。我们对实验进行了 100 次模拟，平均先验转换率为 32%，预期损失阈值𝜖为 0.0015，相对最小可检测效应为 15%(因此我们寻找𝜆_𝑡≥0.15∫𝜆_𝑐).

现在让我们来研究一下模拟数据。

实验模拟(图片由作者提供)

通过检查，我们可以看到，在通过处理达到最小可检测效果的情况下，大多数实验在我们达到每个变体 2000 个样本时结束。然而，我们仍然需要对 4000 个用户进行测试，这是一个很大的数字，可能需要一段时间。因此，让我们看看我们是否可以减少这一点，同时仍然非常肯定，测试将是决定性的。

为了停止测试，我们需要设置一个最小样本数，这样我们就不会掉进偷看的陷阱。让我们放大上面的预期损失图来选择这个最小样本数。

实验模拟—放大低样本量(图片由作者提供)

再次，通过仔细观察，我们可以看到，如果我们将最小样本量设置为 200，大多数假阳性(选择控制的预期损失低于阈值的情况)都可以被消除。我们现在继续，并绘制在每个模拟中需要多少样本来宣布获胜者[2]，假设每个变体至少观察到 200 个样本。

结论性模拟(图片由作者提供)

我们看到，我们只需要每个变体 450 个样本，80%的测试就能得出结论。因此，如果我们决定运行测试，直到我们有 900 个样本，然后，鉴于最低可检测的影响，我们有 80%的机会，测试将是决定性的。因此，我们看到，我们很可能找到一个比我们原来认为的少 3100 个样本的结论性结果。

应用同样的逻辑，现在让我们看看我们前面讨论的差异因素是如何影响所需样本量的。

最小可检测效应

最小可检测效应𝛿是我们想要检测的𝜆的最小相对百分比变化。这是最小的变化，将使治疗值得推出给我们的用户群。让我们看看这个𝛿如何影响我们的样本量计算。当从{0.05，0.15，0.3}中选择𝛿时，我们将使用 32%的先前平均转换率和 0.0015 的预期损失阈值。

𝛿变化的决定性模拟(图片由作者提供)

我们发现，最低可检测效应越高，80%的决定性试验所需的每个变异体的样本量越低，反之亦然。选择 5%的𝛿将要求我们总共获得大约 2100 个样本，而选择 30%的𝛿将要求我们总共只获得大约 600 个样本。直觉上，这是有意义的，因为转换概率的变化越大，我们需要确定它的样本就越少。

预期损失阈值

预期损失阈值𝜖是在我们错误地选择了一个变量的情况下，我们愿意接受的最大预期损失。这是我们在这种情况下最大的预期转化率下降。让我们看看这个𝜖如何影响我们的样本量计算。当从{0.0005，0.0015，0.003}中选择𝜖时，我们将使用 32%的先验平均转化率和 15%的最小可检测效应。

𝜖变化的决定性模拟(图片由作者提供)

我们看到，预期损失阈值越高，80%的检验所需的每个变量的样本量越低，反之亦然。选择 0.0005 的𝜖将需要我们总共获得大约 1350 个样本，而选择 0.003 的𝜖将需要我们总共仅获得大约 650 个样本。所需样本量的变化小于最小可检测效应的情况，因为𝜖的变化规模不同于𝛿.的变化规模然而，概念是相似的，预期损失阈值越低，我们就越想确定实验的结果，因此所需的样本量就越大。

转换概率的标度

现在让我们看看测试中使用的转换概率的比例如何影响我们的样本大小计算。我们将使用 15%的𝛿，并选择 0.005 的相对𝜖。我们在这种情况下使用了相对𝜖，因此分析更加公平，这与在之前的情况下使用𝜖 = 0.0015 相同，之前的平均转换率保持不变，为 32%。在这种情况下，我们将考虑{0.05，0.32，0.6}之前的平均转换率。

𝜆变化的决定性模拟(图片由作者提供)

我们看到一个有趣的结果。尽管对于较高的转换率，所需的样本量没有(相对)大的差异，但是随着转换率变低，所需的样本量变得不成比例地大。在计算所需的样本量时，这是一个需要牢记的重要概念。这是因为转换概率𝜆越低，𝜆的后验分布就越分散，因此我们需要更多的样本来减少这种分散。

我希望这篇文章对你估计贝叶斯产品实验的合理测试持续时间有所帮助，并理解影响它的因素。观看这个空间的系列的下一部分！

参考

【1】贝叶斯 AB 测试对偷看免疫吗？作者大卫·罗宾逊

[2] 贝叶斯 A/B 测试——布莱克·阿诺德的模拟实践探索——我从布莱克的帖子中获得了查看结论性模拟百分位数的想法

我这篇文章的代码可以在这里找到

贝叶斯 AB 检验——第四部分——选择先验

原文：https://towardsdatascience.com/bayesian-ab-testing-part-iv-choosing-a-prior-5a4fe3223bfd?source=collection_archive---------16-----------------------

如何为您的测试度量选择一个合理的先验分布

杰森·登特在 Unsplash 上拍摄的照片

系列结构

这篇文章是关于将贝叶斯 AB 测试方法应用于现实生活产品场景的系列博客文章的第 4 部分。它使用了本系列第一部和第二部中讨论的一些概念。

1. 基于转换的测试指标(速率指标)的建模和分析
2. 基于收入的测试指标(持续指标)的建模和分析
3. 计算测试持续时间
4. 选择合适的先验
5. 用多个变量运行测试

什么是先验分布？

在贝叶斯推理中，先验分布是一种概率分布，用于在从基础总体中抽取样本之前表明我们对未知变量的信念。然后，我们使用这些数据，通过贝叶斯法则更新我们对变量的信念，从而得到变量的后验分布。

在 AB 测试的上下文中，先验分布是一组我们认为测试度量要采用的值，并为每个值分配一个概率。然后，我们以随机实验的形式抽取样本，用于计算后验分布。这些后验分布反过来用于计算 AB 检验的结果。

先验分布的选择如何影响实验？

嗯，贝叶斯法则告诉我们

用文字来说，可以写成

其中分母是归一化常数。因此，该规则可以简化为

因为测试的结果是根据后验概率计算的，而先验是后验概率的一个因素，所以先验的选择对测试有影响，但我们需要小心它的影响有多大。如果我们选择太强的先验，先验将成为主导因素，提取样本的可能性不会产生太大的影响，使实验变得无用。这可能会导致后验概率、控制概率和变异概率快速收敛，从而导致测试不确定。然而，如果我们选择一个非常弱的先验，后验概率将主要依赖于可能性，所以我们需要更多的样本来得出结论性的结果，导致测试时间更长，产品迭代速度更慢。

为了使我们的后验分布更容易计算，我们也可以使用共轭先验。共轭先验是我们可以与似然函数一起使用的先验分布，使得我们计算的后验分布具有与先验分布相似的形式。使用共轭先验简化了我们的计算，同时仍然为测试度量提供了良好的统计模型。在本系列的第一篇和第二篇文章中，我们已经看到了简化的计算和共轭先验的选择是如何对我们有利的。

前科类型

在我们深入探讨如何选择一个先验之前，让我们快速看一下三种主要类型的先验。[1]

主观的

• 基于实验者对该领域的了解
• 在我们的例子中，这将基于产品和数据团队先前对这个测试指标的经验

客观且信息丰富

• 基于该值的历史数据
• 在我们的例子中，这将基于我们所拥有的关于测试指标的任何历史数据
• 它也可能是以前实验的后验分布

无信息的

• 不传达任何价值信息的先验
• 在我们的例子中，这将是测试度量空间上的均匀分布

让我们假设我们是公司和产品的新手，所以我们没有足够的信息来使用主观先验。我们也不希望使用无信息的先验，因为我们认为这将导致更长的测试，从而阻碍我们产品的进展。让我们来看几个技巧，我们可以用来选择一个客观的和信息丰富的先验。

检查

选择先验分布的最简单方法是绘制和检查相关测试指标的历史数据。为了进一步理解这一点，让我们考虑这个系列的第篇文章中的实验。假设我们最近更改了追加销售屏幕上的信息，并希望在向更广泛的用户群发布之前进行 AB 测试。我们假设我们所做的改变将会带来更好的转化率。

在我们设置测试之前，我们希望使用历史数据来选择一个先验。让我们看看如何绘制数据来帮助我们选择。我们将数据分成 100 个分区，计算出每个分区的转换率，并将转换率绘制成直方图。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

prior_data = pd.read_csv('prior_data_conversions.csv')x = np.linspace(0,1,1000)partitions = np.array_split(prior_data, 100)rates = []for partition in partitions:
    rates.append(partition['converted'].mean())_, ax = plt.subplots()sns.histplot(rates, kde=True, label='CR')ax.legend()
ax.set_xlabel('Conversion Rate')
ax.set_ylabel('Density')
ax.set_title('Histogram of Prior Conversion Rates')

先前转换率的直方图(图片由作者提供)

我们现在可以选择一个类似于上面分布的先验分布，但是稍微弱一点。我们不想选择太强的先验，因为我们希望可能性成为计算先验的主导因素。然而，我们确实希望选择一个足够强的先验，这样测试持续时间将会更短。

我们将使用β分布来模拟我们的转换率，因为它是[0，1]上的灵活分布，也是一个很好的共轭先验。所以让我们继续，为我们的练习绘制一些不同强度的潜在先验。

import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt_, ax = plt.subplots(1, 1)x = np.linspace(0,1,1000)beta_weak = beta(4, 8)
beta_mid = beta(16, 36)
beta_strong = beta(33, 69)ax.plot(x, beta_weak.pdf(x), label=f'weak Beta({4}, {8})')
ax.plot(x, beta_mid.pdf(x), label=f'mid Beta({16}, {36})')
ax.plot(x, beta_strong.pdf(x), label=f'strong Beta({33}, {69})')ax.set_xlabel('Conversion Probability')
ax.set_ylabel('Density')
ax.set_title('Choice of Priors')
ax.legend()

先验选择(图片由作者提供)

我们看到，即使我们绘制的最强先验也比转化率的历史分布弱。所以我们可以提前选择𝐵𝑒𝑡𝑎(33,69)作为我们的前期分销。

我们现在可以运行我们的实验，计算后验概率和测试结果。要了解更多关于如何做到这一点，特别是对于概述的实验，请查看这篇帖子。

蒙特卡罗马尔可夫链(MCMC)

一个更复杂但非常有趣的选择先验分布的方法是使用蒙特卡罗马尔可夫链。这种方法对于未知变量由其他随机变量决定的模型特别有用，每个随机变量都有自己的分布。因此，这是一个很好的 AB 测试技术，测试指标是基于收入的(比如每个用户的平均收入)。

在我们开始讨论如何使用这种方法之前，让我先介绍一下它是如何工作的——MCMC 应该有自己的一篇文章，所以这篇介绍会非常简短。[2] MCMC 方法允许我们通过运行模拟从一个未知分布中取样(因此名称的 蒙特卡罗 部分)，其中我们创建一个马尔可夫链，它将我们的未知分布作为其平稳分布。

但是这些术语实际上是什么意思呢？嗯，马尔可夫链是一个在一组状态之间跳跃的过程，每次跳跃都遵循马尔可夫性质。简而言之，这意味着跳转到某个特定状态的概率只取决于流程的当前状态，而不是流程跳转的前一个状态。由于这种无记忆的特性，以及在不同状态之间跳跃的概念，这个过程通常被称为随机行走。让我们假设我们执行无限步的随机行走，那么平稳分布是我们访问每个状态的步数的比例。

现在我们对 MCMC 方法有了一点了解，让我们开始使用它们来为 AB 测试选择先验。让我们考虑一下这个系列的第二篇文章中的实验。我们最近对应用程序中的商店功能进行了 UX 更改。我们相信这些变化让我们的用户更容易进行更大的应用内购买，我们希望在向更广泛的用户群发布之前进行 AB 测试。我们假设我们所做的改变将导致每用户平均收入的显著提高。

我们将每个用户产生的收入建模为随机变量𝑅=𝑋∗𝑌，其中:

𝑋是一个伯努利随机变量，指的是用户是否进行了购买，转换概率为𝜆-𝑋∼𝐵𝑒𝑟(𝜆)

𝑌是一个指数随机变量，指的是购买的规模，如果它与利率参数𝜃-𝑌∼𝐸𝑥𝑝(𝜃)

我们可以用𝜆和𝜃的共轭先验来简化我们的计算。

我们现在需要为我们的参数选择先验，这可能是无信息的。

import arviz as az
import pymc3 as pmprior_revenue = pd.read_csv('prior_data_revenue.csv')rev_observed = prior_revenue[prior_revenue['converted'] == 1]['revenue'].valuesconv_observed = prior_revenue['converted'].valuesmodel = pm.Model()with model:
    alpha = pm.Uniform("alpha", lower=0, upper=100)
    beta = pm.Uniform("beta", lower=0, upper=100)
    k = pm.Uniform("k", lower=0, upper=5)
    theta = pm.Uniform("theta", lower=0, upper=5) cr = pm.Beta('cr', alpha=alpha, beta=beta)
    rr = pm.Gamma('rr', alpha=k, beta=(1/theta)) conversion = pm.Bernoulli('conversion', p=cr, observed=conv_observed) revenue_per_sale = pm.Exponential('revenue_per_sale', lam=rr, observed=rev_observed) trace = pm.sample(10000, return_inferencedata=False)

一旦我们拟合了模型，我们现在可以绘制每个参数的分布，并打印出一些汇总统计数据。

with model:
    az.plot_trace(trace, compact=False)

MCMC 跟踪(图片由作者提供)

with model:
    display(az.summary(trace, kind='stats', round_to=2))

map_estimate = pm.find_MAP(model=model)print(map_estimate)

我们将要使用的两个主要统计数据是每个参数的平均值和每个参数的 MAP 估计值。简而言之，后者是对每个参数分布的点的估计，这些点导致转换率和收益率分布的模式。由于我们的参数先验是一致的，这些估计也是𝜆和𝜃.先验分布的极大似然估计[3]

让我们继续使用这些统计数据来绘制前科。

from scipy.stats import betacr_prior_mean = beta(33, 67)
cr_prior_map = beta(47, 100)x = np.linspace(0,1,1000)_, ax = plt.subplots()sns.lineplot(x=x, y=cr_prior_mean.pdf(x), label='mean Beta(33,67)')
sns.lineplot(x=x, y=cr_prior_map.pdf(x), label='map Beta(47,100)')ax.set_xlabel('Conversion Probability')
ax.set_ylabel('Density')
ax.set_title('Conversion Probability Prior')
ax.legend()

在转换概率𝜆的情况下，两种分布非常相似。我们将继续选择较弱的一个作为度量，因此我们的先验由下式给出

𝜆∼𝐵𝑒𝑡𝑎(33,67)

from scipy.stats import gammarr_prior_mean = gamma(a=2.3, scale=2.0)
rr_prior_map = gamma(a=5, scale=0.4)x = list(range(20))rr_mean = [rr_prior_mean.pdf(i) for i in x]
rr_map = [rr_prior_map.pdf(i) for i in x]_, ax = plt.subplots()sns.lineplot(x=x, y=rr_mean, label='mean Gamma(2.3,2.0)')
sns.lineplot(x=x, y=rr_map, label='map Gamma(5,0.4)')ax.set_xlabel('Revenue Rate')
ax.set_ylabel('Density')
ax.set_title('Revenue Rate Prior')
ax.legend()

类似地，在𝜃收益率的情况下，让我们继续选择较弱的先验，它使用来自我们的 MCMC 算法的𝑘和θ分布的平均值。所以我们有

𝜃∼𝐺𝑎𝑚𝑚𝑎(2.3,2.0)

现在我们有了先验知识，我们可以进行实验，计算后验概率和测试结果。为了找到更多关于如何做到这一点的信息，特别是对于这个概述的实验，请查看这篇帖子。

我希望您发现这种用于选择先验知识的技术探索是有帮助的。请关注本系列的下一部分！

参考

[1]http://www.stats.org.uk/priors/Bayes6.pdf

[2]Rahul Agarwal 的《每个人的 MCMC 直觉》——我发现这对理解 MCMC 算法很有帮助

[3] 最大似然估计 VS 最大后验 A杨

我在这篇文章中的代码可以在这里找到。

通过 JAX 的贝叶斯近似神经网络示例

原文：https://towardsdatascience.com/bayesian-approximated-neural-network-example-via-jax-flexible-parameter-distributions-for-odes-9b4a7f79b493?source=collection_archive---------21-----------------------

常微分方程的灵活参数分布

虽然存在几种已知的方法来近似 ode 的参数分布(特别是 PyMC3 )，但这种概念验证方法在不提供先验分布的情况下估计未知分布。使用神经网络中的漏失概念作为模型不确定性的贝叶斯近似形式，可以通过训练的神经网络的采样来近似灵活的参数分布。

由于具有由给定神经网络学习的灵活输出分布的优势，理论上预测可以更好地拟合伴随有模型不确定性而不是硬离散数字的真实生活数据。对于许多业务场景，尤其是与医疗保健行业的患者护理相关的业务场景，拥有这些预测信息是必不可少的。毕竟，当模型的可解释性或置信度未知时，医生为什么要根据给定患者的图表考虑稳健模型的输出以获得更好的药物治疗呢？在做出关键决策时，我们当然希望对模型的预测有所了解。当将这种方法与常微分方程(ODE)甚至常微分方程系统相结合时，可以更容易更有效地对复杂的真实世界动力学进行建模，而同时仍然应用一些刚性结构(以增加可解释性)。可以导出模型不确定性和灵活的 ODE 参数分布，而不必对潜在参数“可能”是什么做出有影响的先验假设。

启发我的是这篇论文:

【arxiv.org 【1506.02142】作为贝叶斯近似的辍学:表示深度学习中的模型不确定性】

本文中的示例用于估算一个简单 ODE 所需的两个参数，dy/dt = -r * y。一个参数用于增长率 r，另一个参数用于对某个 t 进行积分的初始 y 值。概念验证的一个期望是可视化 t 轴上的输出分布，并有可能在 t 轴上具有可变分布。通过对经过训练的贝叶斯神经网络进行采样(通过丢弃)，我们可以预期近似的 ODE 参数分布是灵活的，因此不能保证是一致的。我发现这对于更复杂的情况是非常强大和有用的，而不必提供一些假设的先验分布，因为神经网络将近似参数分布，并且还可以扩展以及支持连续/在线训练。

我更喜欢将变量依赖关系保存在外部，而不是内嵌，但是我想提供一个工作的单个文件脚本，以便更容易地查看所有引用的函数和变量。下面是单个文件中的完整代码(文件:Bayesian _ appx _ nn _ ex _ jax . py):

https://github.com/bmurders2/Bayesian-Approximated-Neural-Networks-for-Flexible-ODE-Parameter-Distributions

Python 文件依赖于以下包:

• numpy
• 熊猫
• matplotlib
• 海生的
• jax[cpu] ( JAX )

我在执行 Python 脚本时的方法依赖于 Tensorflow GPU Docker 镜像(版本 2.6.0，Python 3.6.9)进行加速的 GPU CUDA 训练。然而，运行这个脚本既不需要 Tensorflow 也不需要 GPU，因为标准的基于 CPU 的 JAX 包也可以工作(用列出的包和 Python 3.9 成功测试)。JAX 对 XLA 和 JIT 的利用甚至对于 CPU 来说都是惊人的！如果您在执行这个 Python 文件时碰巧使用了 VS 代码，那么有一些字符串，#%%，它们是 VS 代码的 Python 交互特性的一部分，用于在处理 Python 文件时获得类似 Jupyter 的体验。如果这个 Python 文件作为 Jupyter 文件导出，这些字符串将自动创建它们相应的单元格(如果直接在 Jupyter 笔记本中运行这个带有交互单元格的文件比 VS Code 的交互窗口特性更好)。

JAX 的实验包包括 Stax 直接在 JAX 境内进行神经网络开发和训练。其他强大的软件包使用 JAX 创建非常健壮的神经网络，如 Trax ，但我认为 Stax 更容易为这个演示建立原型，不需要 JAX 以外的其他软件包。Python 文件使用以下神经网络超参数:

• 辍学率值:0.1
• 每层单位:2048
• 隐藏层数:1
• 主要激活功能:MISH

虽然“深”神经网络将允许更鲁棒的拟合，但是这个示例 ODE 对于单个隐藏层来说足够简单，但是具有“非常宽”的隐藏层(当试图通过神经网络放弃这种贝叶斯近似方法来关联或“近似”高斯过程时，这是一个重要的概念)。对于训练，在对训练数据进行全面训练之前，有一个“好的”猜测是非常重要的。这类似于大多数优化器或最小化算法，其中初始参数值的“坏”猜测可能阻止适当拟合的有效收敛。因此，训练被分成两部分，第一部分用于在对添加了噪声的训练数据进行最终训练之前“启动”神经网络以输出“好的”猜测参数值。采用这种方法将大大有助于在训练期间避免 NaN 值，因为给定函数的参数可能具有非常陡的梯度，这取决于参数所在的位置。为两个训练实例选择的优化器是 JAX 的 ADAMAX，其步长为 2e-3 用于启动训练会话，2e-5 用于针对添加噪声的训练数据的训练。对于损失训练函数，使用 Huber。

在对添加了噪声的训练数据进行训练之后，使用 200 个样本来确定使用任意预测区间的分位数的输出分布。在这种方法中，使用直接标准差+平均值上的分位数绘制预测区间更有用，因为您不必担心区间会忽略 ODE 的重要方面，例如跨越目标 ODE 函数的水平渐近线。 Seaborn 和 Matplotlib 用于以下地块:

• 导出跨 t 轴的输出样本的密度视觉。
• 样品的外观。
• 单个样本的图。
• 近似 ODE 参数的分布。

按作者分类的图像—t 轴上输出分布的密度图

按作者分类的图像-输出的等值线视图

按作者分类的图像-带有预测间隔的单个样本视图

按作者分类的图像-参数 y(0)的样本分布

按作者分类的图像-参数 r 的样本分布

我们可以从近似 ODE 参数的密度图中看到，它们并不完全一致，这正是我在训练后想要看到的。目标 ODE 并不太复杂，但是演示这种方法如何灵活地获得非均匀分布才是目标。重要的是，绘制更多的样本可以影响 ODE 参数估计的可视化近似密度图，但是对于 200 个样本，对于近似 ODE 参数密度“是什么”应该是足够的。

我们可以使用某种形式的增广神经常微分方程，通过这里提出的贝叶斯扭曲，直接逼近该函数。但是，如果我们碰巧对潜在的动态有一些合理的知识或假设，那么利用这种方法会更好。试图用黑盒模型直接模拟潜在动态的一个可能的结果是，我们将无法适当地欣赏参考常微分方程确实存在的水平渐近线。然而，这种方法允许利用神经网络来“学习”目标 ODE 参数的分布，并给出有用的预测区间和中值输出。我对这种方法很感兴趣，并希望其他人会受到启发，利用概率输出进行数据建模(特别是像神经网络这样的黑盒模型)。

PyMC3 中的贝叶斯分层建模

原文：https://towardsdatascience.com/bayesian-hierarchical-modeling-in-pymc3-d113c97f5149?source=collection_archive---------2-----------------------

贝叶斯统计

治愈你模特的健忘症

在 Unsplash 上由 Rob Pumphrey 拍摄的照片

两种方法的故事

有时，当你试图在一个完全异构的数据集上进行机器学习时，你会面临选择的痛苦。作为我所说的异构数据集的一个例子，考虑一个简单的身高数据集。

图片由作者提供。

这个数据真正重要的特征是身高。然而，我们也可以在这里找到两个不同的群体:男性和女性。仅由男性组成的子数据集本身是同质的，仅由女性组成的子数据集也是如此，但是将两者放在一起会产生异质数据集。

给定一些异构数据，您可以:

• 在完整的数据集上构建一个大模型，或者
• 在数据集更小、更同质的部分上构建多个模型。

图片由作者提供。

利弊

构建一个大模型，也被称为**(完全)池化**，通常是最简单的方法:你把所有的样本放在一起，忘记不同的组。然而，如果大模型过于简单，它可能会忽略数据中的不同细微差别，从而导致不符合。当使用高度可解释的模型(如线性模型)时，这个问题会很快发生。如果您可以并且想要使用黑盒方法，如梯度增强，该模型可以自己发现并学习不同的子数据集，但代价是降低了可解释性。

建立几个更小的模型似乎是一个更好的主意:你将数据集分割成不同的更同质的部分，并为每个部分建立一个更小的模型。由此产生的模型被称为非池化模型。这样，每个模型都可以很好地处理一小部分数据，并且它们应该一起形成一个不错的模型。这种方法的明显缺点是你必须适应许多模型，但是当小数据集变得太小时，一个更严重的问题出现了:过度适应。

为了说明这一点，假设我们的身高数据集由 100 个样本组成:97 名女性和 3 名男性。我们现在想通过取身高的样本平均值来估计每组的平均身高。虽然女性的猜测应该是相当不错的，但男性的结果将非常不可靠，因为观察的数量很少。主要问题如下:

子模型是完全独立的。被训练的每个子模型不向其他子模型传递任何信息。在某种意义上，总模型忘记了它在每个子数据集上做什么。

未冷却的模型。

模特的健忘症

虽然训练孤立的模型在某些情况下可能有意义，但在我看来，子模型共享知识往往更好。

举个现实生活中的例子，在 小酒馆到德国金色 M 买个芝士汉堡的等待时间平均需要两分钟(这个数字是我编的)。如果我现在去法国，平均等待时间很可能与这两分钟相差不远。日本可能也是如此。当然，一些国家更快，一些国家平均更慢，但是不管是哪个国家，等待一个汉堡一个小时应该是一个相当罕见的例外。所以，如果我饿了，想在法国吃一个芝士汉堡，我会想到要等多久，也许我是对的。我可以将从一个小组(德国)学到的知识转移到另一个小组(法国)进行预测。这同样适用于 heights 数据集。

一个有用的默认

当然，你也可以创造一些情景，让人们误以为记住了你学过的东西。在医生办公室的等待时间可能不会给你任何关于交通灯绿色阶段持续时间的指示，所以你不应该试图在这两个不相关的变量之间分享知识。

然而，如果我们诚实，通常我们不会看到来自完全不相关和随机来源的数据。通常情况下，这是相当一致的，只是在国家、性别、商店等群体中有一些变化，因此

我觉得分享知识往往是一个很好的默认。

那么，如何让我们的模型有更好的记忆呢？或者说，换一种方式，如何教会子模型之间的沟通和协同工作？这有几种方法，例如，在神经网络的情况下，参数或权重分配，这里我不打算详细说明。

然而，正如你从这篇文章的标题中所知道的，我们将转向贝叶斯，并学习如何做贝叶斯分层建模作为我们问题的可能解决方案。我们将再次使用 PyMC3 ，因为它是一个不错的包。

如果您从未听说过贝叶斯统计或 PyMC3，并想学习它，请查看我关于这个有趣主题的其他介绍性文章。

走吧！

PyMC3 中的分层建模

首先，我们将重新审视两者，在贝叶斯设置中的池化和非池化方法，因为它是

1. 一个很好的练习，而且
2. 非池化和层次化(也称为部分池化或多级)的代码库非常相似。

在我们开始之前，让我们创建一个数据集进行实验。我们将创建一个简单的一维回归问题，即只有一个特征和一个目标。有八个不同的组，每个组都有自己的斜率，固定截距为零。

import numpy as np

np.random.seed(0) # to keep it reproducible

mean_slope = 2 # the 8 different slopes have a mean of 2

slopes = np.random.normal(mean_slope, size=8)
groups = np.array(50*[0, 1, 2, 3, 4, 5, 6] + 5*[7])

x = np.random.randn(355)
y = slopes[groups] * x + 0.1*np.random.randn(355)

groups变量包含每个观察值所属的组。0 到 6 组各有 50 个观察值。然后，最后是少数群体 7，只有五个观察值。

作为一个小指南，要获得属于组 2 的元素，可以使用x[groups==2], y[groups==2]。这个子数据集的斜率是slopes[2]。

为了增加趣味，让我们在这个少数群体中加入局外人。

y[-1] = 30 # the last 5 observations are from minority group 7
y[-2] = 30

导入我们最喜欢的贝叶斯库后

import pymc3 as pm
import arviz as az

我们可以开始建模了！

混合模型

我们将从忽略组开始，将这个数据集视为一个大块。PyMC3 中不带截距的简单贝叶斯线性回归如下所示:

with pm.Model() as pooled_model:
    slope = pm.Normal('slope', 0, 20)
    noise = pm.Exponential('noise', 0.1)

    obs = pm.Normal('obs', slope*x, noise, observed=y)

    pooled_trace = pm.sample(return_inferencedata=True)

az.plot_posterior(pooled_trace, var_names=['slope'])

输出将是

图片由作者提供。

这基本上没用，因为我们只得到一个单一的斜率，但这并不奇怪。下一个最好的方法是为每组引入一个斜率。

无池模型

with pm.Model() as unpooled_model:
    slope = pm.Normal('slope', 0, 20, shape=8)
    noise = pm.Exponential('noise', 10)

    obs = pm.Normal('obs', slope[groups]*x, noise, observed=y)

    unpooled_trace = pm.sample(return_inferencedata=True)

az.plot_posterior(unpooled_trace, var_names=['slope'])

我们现在有八个不同的斜坡，都是独立训练的:

图片由作者提供。

我们也可以仅标绘 94%的高密度区间(HDI ),即包含 94%后位质量的短可信区间

az.plot_forest(unpooled_trace, var_names=['slope'], combined=True)

我们得到了

图片由作者提供。

你可以看到第 0 组到第 6 组的斜率很小，而第 7 组的斜率很大。但这是完全错误的，因为我们所有的斜率都应该在值2左右。发生了什么事？简单:我们通过在最小的群体中引入异常值来欺骗模型。

这正是我之前谈到的问题:组 7 的子模型没有机会，因为它不知道组 0 到组 6 中发生了什么。它不知道斜率通常在 2 左右。

所以，让我们来解决这个问题。

部分池化 aka 层次模型

with pm.Model() as hierarchical_model:
    mu_slope = pm.Normal('mu_slope', 0, 1) # hyperprior 1
    sigma_slope = pm.Exponential('sigma_slope', 13) # hyperprior 2

    slope = pm.Normal('slope', mu_slope, sigma_slope, shape=8)
    noise = pm.Exponential('noise', 10)

    obs = pm.Normal('obs', slope[groups]*x, noise, observed=y)

    hierarchical_trace = pm.sample(
        return_inferencedata=True,
        target_accept=0.995
    )

az.plot_posterior(hierarchical_trace)

那么，现在这是什么？在无池模型中，我们通过

slope = pm.Normal('slope', 0, 20, shape=8)

我们告诉模型，斜率应该在零左右，但是具有相当大的 20 的标准偏差。

在分层模型中，我们定义了所谓的超先验来为斜率先验找到更好的均值和标准差。重要的是:

mu_slope = pm.Normal('mu_slope', 0, 1)
sigma_slope = pm.Exponential('sigma_slope', 13)

slope = pm.Normal('slope', mu_slope, sigma_slope, shape=8)

我们用mu_slope代替零，用sigma_slope代替二十，就这么简单。这两个都是随机变量，我们可以用贝叶斯推理来学习。mu_slope和sigma_slope被称为超优先级，例如，决策树的最大深度被称为超参数。它们都比slope高一个层级，因为在slope可以计算之前，它们必须首先被评估。

注意:我在sample方法中添加了target_accept=0.995来改进采样，因为从这个嵌套的后验模型中采样不再像对非 pool 模型那样容易。我可以在另一篇文章中进一步解释这一点。

后验看起来像这样:

图片由作者提供。

图片由作者提供。

这仍然不是完美的，但比无池模型更接近事实。这是因为该模型还试图使用mu_slope超级函数来计算斜率平均值。真实的先验应该是 2，然而，由于异常值，模型认为它大约是 3。

图片由作者提供。

尽管如此，我们在非 pool 模型中看到的大斜率15将拉向这个大约 3 的mu_slope。其他参数也是如此:

在分层建模中，所有参数都被拉向全局平均值。这种效应被称为收缩。

总而言之，层次模型是明显的赢家。作为最后一步，让我们使用标牌符号直观地比较三种型号。这是一个不言自明的模型可视化表示，可能对你们当中的视觉学习者有所帮助。

使用平板符号比较模型

PyMC3 有一个奇妙的功能pm.model_to_graphviz(model)以一种很好的方式显示你的模型。

图片由作者提供。

我喜欢这种表示，因为你可以清楚地看到模型的某种演变。在混合模型中，我们从单一斜率开始。然后我们增加斜率的数量，用数字为 8 的方框表示。

这个符号是为斜坡画 8 个圆的捷径。我猜这个小盒子应该是一张桌子，上面叠放着 8 个盘子，我们从上面看，因此得名。

这个过程的最后一步是在顶层引入一个新的变量层，即超优先级。

结论

在本文中，我们讨论了处理异构数据集时不同的建模方法。起初，我声称有两种方式:池化和非池化模型。很明显，我在撒谎，正如我们已经看到的分层建模一样。

一个有趣的观察是，您可以将分层建模看作是对非池化建模的推广。如果您将超优先级设置为某个常量随机变量，您将再次以无池方法结束。否则，层次模型通过跟踪全局参数来做一些不同的事情，例如mu_slope，它充当我们想要估计的实际参数的种子。作为一个很好的副作用，我们也得到这些全局参数的估计。如果在不同的组中没有值得分享的东西，层次模型也会学习这个，这很好。这就是为什么我认为你应该尽可能使用层次化的方法，而不是非池化的方法**:通常这没有坏处。**

当使用分层建模而不是无池建模时，唯一可能出错的是有更多的参数要估计，并且由于嵌套的参数，估计也变得更加困难。你必须使用一些技巧，比如改变 MCMC 采样参数或者引入所谓的非中心变量。

一个以为中心的 ( 我称之为嵌套的)变量类似于另外两个随机变量a和b的pm.Normal('X', a, b)，这就是我们上面用过的。然而，在正态分布的特殊情况下，你也可以写a + b*pm.Normal('X', 0, 1)，这从统计学的角度来看是等价的，但对于计算后验分布的 MCMC 算法来说是一个巨大的差异。

尽管如此，如果你能克服计算问题，我认为分层建模通常是一条值得尝试的道路。

作为一个额外的资源，看看托马斯·威奇和达内·埃尔伯斯写的这篇不错的博文。

我希望你今天学到了新的、有趣的、有用的东西。感谢阅读！

作为最后一点，如果你

1. 想支持我多写点机器学习和
2. 无论如何，计划获得一个中等订阅量，

为什么不做 通过这个环节 ？这将对我帮助很大！😊

透明地说，给你的价格不变，但大约一半的订阅费直接归我。

非常感谢，如果你考虑支持我！

如果你有任何问题，请在 LinkedIn 上给我写信！

Python 中的贝叶斯推理和马尔可夫链蒙特卡罗抽样

原文：https://towardsdatascience.com/bayesian-inference-and-markov-chain-monte-carlo-sampling-in-python-bada1beabca7?source=collection_archive---------3-----------------------

通过一个用 Python 实现的掷硬币的深入示例，介绍如何使用贝叶斯推理和 MCMC 采样方法来预测未知参数的分布。

图片来自 Adobe Stock

介绍

本文将一个基本的抛硬币例子推广到一个更大的环境中，在这个环境中，我们可以检查贝叶斯推理和马尔可夫链蒙特卡罗抽样在预测未知值方面的用途和能力。有许多有用的软件包可以使用 MCMC 方法，但是这里我们将从头开始用 Python 构建我们自己的 MCMC，目标是理解其核心过程。

什么是贝叶斯推理？

贝叶斯推断是一种方法，当我们收集更多的数据和证据时，我们使用贝叶斯定理来更新我们对概率或参数的理解。

什么是马尔可夫链蒙特卡罗抽样？

MCMC 方法(通常被称为)是一种用于从概率分布中取样的算法。这类算法采用随机抽样来获得数值结果，随着样本数量的增加，这些结果会收敛于事实。有许多不同的算法可以用来创建这种类型的采样链——我们在这个具体例子中使用的算法称为 Metropolis-Hastings 算法。MH 实现允许我们从一个未知的分布中抽取样本，只要我们知道一个与我们想要从中抽取样本的分布的概率密度成比例的函数。不需要知道精确的概率密度，仅仅是它的比例性就使得 MCMC-MH 特别有用。我们将在这个例子的后面描述算法的详细步骤。

硬币工厂问题

这个例子将模拟掷硬币的结果的概率，但是在一个扩展的上下文中。想象一下，一个全新的硬币工厂刚刚建成，生产一种新的硬币，他们让你确定工厂生产的硬币的一个参数，即硬币翻转时正面/反面落地的概率。我们称之为硬币的“偏向”。每枚硬币都有自己的偏差，但鉴于它们是在相同的过程中生产的，只有轻微的差异，我们预计工厂不会随机生产硬币，而是围绕“工厂偏差”生产。我们将使用来自特定硬币和工厂偏差的信息来创建概率分布，以确定硬币最有可能产生的偏差。用贝叶斯的话来说，我们将使用可能性(硬币数据)和先验(工厂偏差)来产生硬币偏差的后验分布。

图片来自 Adobe Stock

于是，工厂开张并生产了第一枚硬币。鉴于它是一种新的硬币和铸币厂，我们不知道它在翻转时会有什么表现(假装我们不知道硬币通常是如何工作的，我们试图预测一个“未知”的参数)。因此，我们用第一枚硬币做了一个实验，将它翻转 100 次，记录下它正面着地的次数——57 次。为了更新我们对工厂偏差的理解，我们将在生产更多硬币时这样做，但首先，让我们只分析第一枚硬币。

贝叶斯背景下的问题

在我们走得太远之前，我们需要在贝叶斯定理的背景下定义这个问题。等式的左边称为后验；一般来说，就是一个假设( H )给定一些证据( E )的概率。在右边的分子中，我们有我们的可能性(假设我们的假设为真，看到证据的概率)，乘以先验(假设的概率)。左边的分母是边际可能性或证据；这是观察证据的概率。幸运的是，我们将不需要使用边际可能性来抽样后验概率。

贝叶斯定理-作者图片。

在大多数实际应用中，直接计算后验分布是不可能的，这就是为什么我们采用数字采样技术，如 MCMC。

在后面的

那么我们对硬币厂感兴趣的后验概率是多少呢？这是概率， P ，工厂生产一个有偏差的硬币， p ，给定我们的数据， x 人头数— P(p|x) 。定理的其余部分可以写成如下:

我们硬币工厂背景下的贝叶斯法则。

这里需要记住的是，我们预测的是后验概率分布( P(p|x) )一枚硬币有一定概率正面朝上，或者偏向( p )，这是两回事。工厂偏差是以一定偏差生产的硬币的概率分布；这是 P§ ，先验。

可能性—二项式分布

这里的似然函数是观察到正面的概率， *x，*给定一个有偏差的硬币 p. 对于抛硬币，这个函数可以用二项分布来精确描述。对于有偏差 p 的硬币，观察到 x 正面出 n 翻转的概率可以写成:

二项式分布-作者图片

值得注意的是，我们目前不知道给定硬币的值 p 。这个值是我们的 MCMC 将随机抽样的值。用一个随机值初始化 p 后，经过多次采样，模型会收敛到 p 的真实值。

让我们开始用 Python 编写这个例子。我们需要定义第一次硬币实验中观察到的数据，以及根据给定的数据从二项分布中得出概率的似然函数。Scipy.stats 是一个很棒的 Python 库，可以很容易地定义像二项式这样的分布，这也是我在这个例子中使用的。

在先的；在前的

接下来，我们需要定义我们的先验函数， P§ 。p 的值只能在 0 和 1 之间，0 代表硬币永远不会正面朝上，1 代表硬币永远正面朝上。请记住，工厂只生产了一枚硬币，所以我们没有任何信息可以预期 p 的先验概率是多少(假装我们不知道硬币是如何工作的)。在这个问题的背景下，由于只生产了一枚硬币，我们还不知道工厂偏见可能是什么。因此，我们将使用所谓的均匀先验。在贝叶斯推理的上下文中，这意味着我们赋予概率 p 在 0 和 1 之间的任意值相等的权重。

我们可以用 Scipy 和均匀分布 PDF 轻松做到这一点。这种默认分布只存在于 0 到 1 之间，正如我们所需要的。

既然我们已经定义了似然和先验，我们可以继续理解和编码马尔可夫链。

大都会-黑斯廷斯 MCMC

如上所述，这些方法从连续的随机变量中抽取样本——在我们的例子中是针对 p 。我们将使用的 MCMC 是随机游走类型，它随机生成样本，并根据它们与模型的拟合程度来决定是否保留它们。

接受比率

Metropolis-Hastings 算法相当简单，但首先，我们需要定义如何接受或拒绝新的样本抽取。每次迭代，将为 0 和 1 之间的 p 建议一个新值，我们将这个建议值称为p′。如果这个值比前一个值更好，我们只想接受并更新它。这是通过计算接受率来完成的。这个接受率是我们的贝叶斯定理对建议值与先前值的比率，如下所示。

接受率， R —作者图片

这里有几点需要注意。首先，你可能已经注意到，这个接受率不包括贝叶斯定理的*边际似然(证据)*部分，我们也没有在上面为它定义一个函数。这是因为证据不会因为新的 p 值而改变，因此在这个比率中抵消了它自己。这太棒了，因为计算贝叶斯定理的边际似然部分在实践中通常是极其困难或不可能的。MCMC 和贝叶斯推理允许我们在不知道边际可能性的情况下对后验样本进行采样！

第二，这里任何大于 1 的值都意味着建议值更好，应该接受。接受新值的第二部分是将 R 与 0 和 1 之间的另一个随机抽取值进行比较，因此习惯上只在 R 更高时将 R 设为 1。

为了简单起见，我们将编写一个函数来计算这个接受率，以便在我们的采样链循环中轻松实现。

让我们明确定义 MH 算法的步骤，这样就更清楚了。

采样算法

我们已经定义了可能性、先验和接受概率的函数。在循环之前，我们必须做的最后一件事是用 0 到 1 范围内的随机值初始化 p 。

以下是我们的 Metropolis-Hastings 算法的步骤:

1. 在 0 和 1 之间随机提出一个 p 的新值，称之为p′(或 p_new)。
2. 计算接受率 r。
3. 生成另一个介于 0 和 1 之间的均匀随机数，称之为 u 。
4. 如果 u < R ，接受新值并设置p = p′。否则，保持 p 的当前值。
5. 记录该样品的 p 的最终值。
6. 多次重复步骤 1 到 5。

非常简单。在我们编写代码之前，让我们先来研究一下 MCMC 中使用的其他一些常见概念。

老化

MCMCs 是随机初始化的，必须收敛到正确的值，这通常需要相当多的样本。当绘制我们的结果和后验分布时，在模型收敛之前包含这些早期样本是无效的。因此，我们实施了所谓的“老化”,即排除第一批不太准确的样本。当整个链约为 10k–20k 时，MCMCs 的老化通常约为 2000–5000 个样本。

落后

MCMCs 需要考虑的另一个非常重要的问题是样本独立性。这里的新样本通常依赖于前一个样本，因为有时我们不接受新的随机值而保留旧值。为了解决这个问题，我们实现了所谓的“滞后”。滞后是指我们不是记录每一个样本，而是每隔一个，或者每五个或十个样本记录一次。

模拟

太好了，我们现在拥有了编写和运行 MCMC 所需的一切。

可视化结果

MCMC 结果通常以两种方式绘制——后验分布和迹线图。后验概率可以显示在直方图中，在直方图中我们可以直观地检查最可能的值和方差。轨迹图显示了每个样本迭代的 p 的值，并显示了 MCMC 的行为和收敛。后验分布不应包含老化样本，但包含迹线老化可以帮助我们检查模型从哪里开始，以及它收敛得如何。老化样本不包括在下面的跟踪图中。

后验分布和迹线——作者图片。

检查结果，我们可以看到我们的后验正态分布。我们还可以从跟踪图中看到，我们在收敛值附近进行了很好的随机采样，这很好。当提取后验值时，对于分布的其余部分，使用 p 的最后值并不总是准确的。由于这个原因，后验值通常被认为是分布的平均值。在这种情况下，我们的平均后验值是 0.569。

这个值几乎是我们预测的，如果我们把这个硬币的数据的频率概率，57/100 头。我们的 MCMC 预测这一点的原因是因为这是第一枚硬币，我们不太了解它们应该如何表现。在均匀先验的情况下，后验受似然函数(即数据)的影响更大。贝叶斯推断的下一步是用更多的数据更新我们的理解，所以让我们保持工厂运转，制造更多的硬币来测试。

照片由 GSJJ 在 Unsplash 上拍摄

更新我们的理解

我们等一会儿，工厂生产了 500 枚硬币。我们运行同样的 100 次翻转实验，并记录每枚硬币的最可能偏差。让我们绘制一个直方图，并检查产生的偏差的分布，以了解工厂偏差。

500 硬币工厂偏见-作者图片。

检查这些数据，我们了解到平均偏差是 0.511，标准偏差是 0.049。数据看起来是正态的，所以让我们用这些参数的正态分布对其建模，这显示为红色。

这个分布包含了我们预计在这个工厂生产的硬币更有可能出现哪些偏差的信息。像这样更新我们的理解，将会给我们一个比我们以前的，不知情的，一致的先验更准确的结果。这正是贝叶斯推理的目的——我们可以简单地更新我们的先验函数来表示工厂偏差的数据。

工厂偏见先验-作者图片。

现在我们已经用更多的先验信息改进了我们的模型，让我们生产第 501 枚硬币并进行同样的实验。在这个实验中，我们在 100 次投掷中得到 63 次投掷。我们必须像以前一样根据这些数据建立我们的似然函数:

很好，现在让我们看看包含工厂偏差的信息如何影响 MCMC 的结果，以及我们认为这枚硬币的真实偏差是多少。我已经运行了一个 MCMC，使用了与上面相同数量的样本和参数，但是使用了更新的先前和新的硬币数据。下图显示了第 501 枚硬币的后验分布和迹线。

后验和迹线来自 MCMC，更新了前验——由作者更新。

尽管这枚硬币的数据表明偏差约为 0.63，但我们的贝叶斯推理模型表明实际值更接近 0.53。这是因为我们的知情先验函数在模型中占有权重，告诉我们即使我们观察到这种硬币有 63 个头像，鉴于硬币的平均偏差约为 0.51，我们预计第 501 个头像的偏差更接近工厂偏差。即使这枚硬币有 0.5 的偏差，观察 100 次投掷中的 63 次投掷也不是完全不可能的，我们不应该假设这个数据代表确切的值。先验函数具有权重，就像似然性具有通知后验分布的权重一样。如果我们要生产成千上万的硬币，并告知先前的分布更多，这将在模型中给予它更高的权重。这种用更多信息来更新我们的理解以预测未知参数的想法正是贝叶斯推理有用的原因。正是利用更多更好的数据来调整和操作这些可能性和先验函数，使我们能够改进和充实我们的推理模型。

结论

当实现贝叶斯推理和 MCMC 时，人们通常不希望知道后验分布是什么样子，因为我们在这里是编造数据的。这是推断未知参数的强大工具——只需知道后验分布与什么成比例(例如，线性模型和预测参数斜率和截距，或逻辑模型和参数α和β)。这允许我们做一些事情，比如预测我们无法解决的高维积分的结果——信不信由你，你可以只用随机性和统计学来做！

参考

除了上面给出的维基百科链接之外，这篇文章对于结合贝叶斯推理和 MCMC 的上下文非常有帮助。

接下来，我建议利用 PyMC3 这样的包来进行贝叶斯推理，而不是从头开始编码。Will Koehrsen 的这篇文章提供了一个很棒的真实例子，值得一看:

如果你希望深入研究使贝叶斯推理和 MCMC 成为可能的数学和推理，我强烈推荐这篇文章— 贝叶斯推理问题，MCMC 和变分推理。

如果你想了解更多关于我的实施的信息，或者想给我反馈(非常感谢)，请给我发信息或电子邮件到 jonny.hofmeister@gmail.com。

谢谢大家！

贝叶斯线性回归:基于 arm 的汽车销售分析

原文：https://towardsdatascience.com/bayesian-linear-regression-analysis-of-car-sales-with-arm-in-r-b1147e79a9a4?source=collection_archive---------22-----------------------

使用贝叶斯线性回归解释不确定性

来源:图片来自 Pixabay

线性回归是最常用也是最有用的建模工具之一。

虽然没有任何形式的回归分析可以接近现实，但它可以很好地预测因变量，并确定每个自变量影响因变量的程度，即每个系数的大小和重要性。

然而，传统的线性回归可能有缺点，因为这种方法不能真正解释估算中的不确定性。

然而，贝叶斯线性回归可以作为这个问题的解决方案-通过重复模拟提供许多不同的系数值估计。例如，如果我们希望进行 20，000 次模拟，这种方法可以返回 20，000 条不同的回归线，这些回归线将说明在各种情况下因变量的潜在波动。

为了说明这一点，让我们来看看 R 中的 arm 库如何使用贝叶斯线性回归来预测汽车销量。

标准线性回归

该分析是在 Kaggle 可用的汽车销售数据集上进行的，该数据集根据 CC0 1.0 通用(CC0 1.0)公共领域专用许可提供。

在第一种情况下，设计线性回归来确定销售额(以千计)如何根据以下因素波动:

• 发动机尺寸
• 车型
• 马力
• 轴距
• 宽度
• 长度
• 整备重量
• 燃料容量
• 燃料效率
• 功率性能因数

这里是线性回归和输出结果。

Call:
lm(formula = Sales_in_thousands ~ Engine_size + Vehicle_type + 
    Horsepower + Wheelbase + Width + Length + Curb_weight + Fuel_capacity + 
    Fuel_efficiency + Power_perf_factor)Residuals:
    Min      1Q  Median      3Q     Max 
-104.47  -29.09   -6.86   14.00  337.19Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)           -302.7471   154.7623  -1.956  0.05242 . 
Engine_size             31.7419    11.0641   2.869  0.00475 **
Vehicle_typePassenger  -45.1881    19.3850  -2.331  0.02117 * 
Horsepower              -0.5830     0.8825  -0.661  0.50990   
Wheelbase                3.7231     1.3670   2.724  0.00727 **
Width                    0.8454     2.4341   0.347  0.72889   
Length                   0.7609     0.8295   0.917  0.36054   
Curb_weight            -63.7597    20.4988  -3.110  0.00226 **
Fuel_capacity           -1.1049     2.7603  -0.400  0.68955   
Fuel_efficiency          0.5560     2.3522   0.236  0.81350   
Power_perf_factor        0.2592     1.9605   0.132  0.89502   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 58 on 141 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.3391, Adjusted R-squared:  0.2922 
F-statistic: 7.235 on 10 and 141 DF,  p-value: 3.643e-09

我们可以看到，发动机尺寸、车型、轴距和整备质量在 5%的水平上具有统计显著性(一些变量在 1%的水平上具有显著性)。

在删除无关紧要的变量并在迭代中再次生成回归后，最终的回归结果如下:

Call:
lm(formula = Sales_in_thousands ~ Vehicle_type + Wheelbase + 
    Curb_weight)Residuals:
    Min      1Q  Median      3Q     Max 
-107.15  -32.26   -9.68   18.33  339.38Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -288.840     76.741  -3.764 0.000239 ***
Vehicle_typePassenger  -44.626     12.673  -3.521 0.000568 ***
Wheelbase                5.287      0.816   6.479 1.23e-09 ***
Curb_weight            -57.100     10.697  -5.338 3.39e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 58.67 on 151 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.2787, Adjusted R-squared:  0.2644 
F-statistic: 19.45 on 3 and 151 DF,  p-value: 1.027e-10

重要变量是车辆类型、轴距和整备质量。虽然 27.87% 的 R 平方在统计学上较低，但这并不一定意味着模型不好。事实上，有许多不同的因素可以影响汽车销售，这里没有考虑。

模型中的其余变量被表示为非常重要，这表明它们与确定汽车销售相关。让我们更详细地分析模型系数。

• 在这种情况下， -288.84 的截距是虚假的，但理论上代表了所有其他系数为零时的最小销售额。但是，销售额不能为负，因此系数本身没有任何意义。
• 由于 Vehicle_typePassenger 变量的系数为-44.626，这意味着如果一种车型被归类为“客车”而不是“轿车”，销售额将下降 $44.626 。
• 轴距每增加一个单位，汽车销量就会增加5.287 美元。
• 整备质量每增加一个单位，汽车销量就会减少**-57.10 美元**。

使用方差膨胀因子检验时，没有一个变量的多重共线性检验为阳性(即 VIF 统计值大于 5)。这表明没有一个变量彼此显著相关，因此扭曲了回归结果。

> library(car)
> vif(M1)
Vehicle_type    Wheelbase  Curb_weight 
    1.384649     1.746395     2.035252

直方图表明汽车销售的分布是非正态的。因此，在这种情况下，异方差性没有被测试(或者实际上被校正)。事实上，正态性并不是生成线性回归模型的必要条件。

来源:RStudio 输出

贝叶斯线性回归

“所有的模型都是错的，但有些是有用的。”

乔治·e·p·Box

虽然上述线性回归提供了信息，但它的一个缺点是没有考虑到估计中的不确定性。

例如，我们真的可以期望轴距每增加一个单位，汽车销量总会增加 5.287 美元吗？

事实上，我们知道情况并非如此。相反，我们想想出一种方法来解释估计中的不确定性，并生成一系列广泛的销售预测。

毕竟，贝叶斯分析的目的是以概率的方式分析数据——考虑尽可能多的潜在情况。

考虑到这一点，R 中的 arm 包用于生成 20，000 条不同的回归线，以说明估计中的不确定性。

我们可以看到，生成模拟时，会生成具有不同系数估计的回归。

> library(arm)
> M1.sim <- sim(M1, n.sims=20000)
> M1.sim
> coef.M1.sim <- coef(M1.sim)
> coef.M1.simAn object of class "sim"
Slot "coef":
         (Intercept) Vehicle_typePassenger Wheelbase Curb_weight
    [1,]  -142.29350           -40.1585711  4.046482   -62.73836
    [2,]  -311.43351           -46.1396357  6.056344   -74.03777
...
  [249,]  -249.44178           -52.5240319  5.590081   -76.08825
  [250,]  -184.34866           -60.6832116  4.493316   -59.79492
 [ reached getOption("max.print") -- omitted 19750 rows ]

此外，还会产生残差(由σhat提取):

> sigma.M1.sim <- sigma.hat(M1.sim)
> sigma.M1.sim
[1] 53.53439 56.59501 62.36744 48.25297 62.11458 65.81519 51.64420
[8] 62.34890 55.86902 52.15755 60.12543 59.68573 55.92115 62.00584
[15] 60.30506 63.08531 56.79228 57.74138 53.94584 57.22139 60.00236
...
[981] 54.00598 63.18356 59.70908 60.27791 64.86504 53.71799 59.88211
[988] 58.08702 56.12842 62.05933 61.02171 54.30973 60.43517 61.74605
[995] 54.91794 60.25361 56.72824 57.06674 61.47696 62.64697
[ reached getOption("max.print") -- omitted 19000 entries ]

现在，可以在分位数的基础上计算系数的不确定性:

> apply(coef(M1.sim), 2, quantile)
     (Intercept) Vehicle_typePassenger Wheelbase Curb_weight
0%    -589.44627             -98.60129  1.929788   -97.18706
25%   -341.31489             -53.40231  4.737972   -64.23857
50%   -289.01086             -44.73874  5.287846   -57.07201
75%   -237.68266             -36.07124  5.836498   -49.90656
100%    20.27591              11.40947  8.547013   -19.46040

通过观察这些系数，我们可以看到系数本身可能存在显著的变化。

例如，虽然乘用车类型降低了 0%、25%、50%和 75%分位数的销售额，但在 100%分位数的模拟中，销售额实际上是增加的。有没有逆势而上的特殊车型？这个模型可以为进一步研究这个问题提供一个有用的起点。

对于汽车的轴距(前轮和后轮之间的距离)，我们可以看到，在 0%分位数时，每增加一个单位对销售额的影响是 1.92 美元，但在 100%分位数时会增加到 8.54 美元。

除了系数估计值之外，我们还可以绘制残差图(由 sigma.hat 提取):

来源:RStudio 输出

以下是基于分位数报告的残差:

> q
      0%      25%      50%      75%     100% 
47.52279 56.60056 58.82638 61.15704 78.96791

鉴于我们已经看到汽车销售的分布具有强烈的正偏态分布(即，选定的几款汽车显示的销售数字比其余的高得多)，这可能是我们看到高不确定性估计值的原因之一，即，50%分位数的不确定性估计值为 58，826 美元。

然而，贝叶斯线性回归在通过生成大范围的模拟和潜在情景来模拟这种不确定性时非常有用。

结论

在本文中，您已经看到:

• 传统线性回归的缺点
• 如何生成贝叶斯线性回归
• 使用 arm 包生成许多回归模拟
• 模型估计的不确定性分析

非常感谢您的宝贵时间，非常感谢您的任何问题或反馈。你可以在michael-grogan.com找到更多我的数据科学内容。

免责声明:本文是在“原样”的基础上编写的，没有担保。它旨在提供数据科学概念的概述，不应被解释为专业建议。本文中的发现和解释是作者的发现和解释，不被本文中提到的任何第三方认可或隶属于任何第三方。作者与本文提及的任何第三方无任何关系。

Bambi 贝叶斯线性回归

原文：https://towardsdatascience.com/bayesian-linear-regression-with-bambi-a5e6570f167b?source=collection_archive---------4-----------------------

利用贝叶斯推理得到你的预测分布

动机

用回归线拟合样本数据时，可能会得到如下的回归线:

作者图片

如果你能得到一个预测的分布，而不是一条回归线，不是很好吗？

作者图片

这就是贝叶斯线性回归派上用场的时候。在本文中，您将学习如何用 Python 构建贝叶斯线性回归。

线性回归到底是什么？

要了解贝叶斯线性回归和普通最小二乘线性回归(OLS)(我们都很熟悉的线性回归)的区别，我们先来了解一下什么是 OLS。

从一个例子开始

理解算法最简单的方法是从一个例子开始。

我们将在 Kaggle 上使用显示学生数学和葡萄牙语成绩的数据集。

作者图片

绘制数据:

作者图片

我们的目标是找到最能描述数学成绩和葡萄牙语成绩之间关系的线性回归线。

如果您有𝑛变量，您的线性回归方程将如下所示:

作者图片

其中𝑦是响应变量，𝑥’s 是解释变量，𝛽ᵢ是变量𝑥ᵢ的权重， ε 是误差项。

由于我们的数据中只有一个解释变量，我们的线性回归方程将如下所示:

作者图片

其中，𝑦是预测的葡萄牙语分数，𝑥是观察到的数学分数。

在矩阵形式中，等式变为:

作者图片

𝛽有几个可能的值，但我们将只选择最小化残差平方和(RSS)的𝛽值。RSS 是𝑦ᵢ和𝑦̂ᵢ之间的平方差之和，其中𝑦ᵢ是数据的第一个实际值，𝑦̂ᵢ是使用回归的第一个预测值。

作者图片

最小化 RSS 将给出𝛽.的估计让我们尝试使用 scikit-learn 的线性回归模型来应用这个方法。

绘制回归线:

作者图片

求系数:

array([0.4335796])

找到截距:

7.403346211071996

通过看系数和曲线图，我们可以说，如果数学成绩增加 1，葡萄牙语平均成绩增加 0.43。如果数学成绩为 0，葡萄牙语成绩的预测值为 7.403。

普通最小二乘法返回单点估计。但是如果我们想看到可能的估计值的分布呢？让我们试着用斑比来做这件事。

使用 Bambi 构建贝叶斯线性回归模型

斑比是什么？

Bambi 是用 Python 编写的高级贝叶斯建模接口。它构建在编程框架 PyMC3 之上。因此，Bambi 类似于 PyMC3，但是更容易使用。

要安装 Bambi，请键入:

pip install bambi

建立线性回归模型

使用 Bambi 构建线性回归模型非常简单。既然我们要用mat(数学成绩)来预测por(葡萄牙成绩)，我们可以写:

使用 PyMC3 从 4 条链中各抽取 1000 个样本(单次运行 MCMC ):

默认情况下，Bambi 使用family="gaussian"，这意味着具有正常误差的线性回归。

分析系数的分布

使用 ArviZ 绘制结果:

作者图片

在上图中，Intercept和mat分别是来自方程𝑦 = 𝛽₀ + 𝛽₁𝑥 + *ε的𝛽₀和𝛽₁。*左边的图显示了变量的分布。右边的图显示了每次迭代中每个变量的值。

看左边的图，我们可以说:

• 𝛽₀和𝛽₁的分布接近正态分布。
• 𝛽₀的平均值从 7.2 到 7.5。
• 𝛽₁的平均值从 0.42 到 0.44。

这太酷了！我们得到的不是𝛽₀和𝛽₁的单点估计，而是𝛽₀和𝛽₁的分布。

分布的形状也有助于我们理解我们对𝛽₀或𝛽₁.的价值有多确定具体来说:

• 分布的扩散越大(越大*【σ】*)，我们对𝛽₀和𝛽₁.的价值就越不确定
• 分布的扩散越小(越小*【σ】*)，我们对𝛽₀和𝛽₁.的价值就越有把握

作者图片

比较贝叶斯线性回归和 OLS 的结果

与 scikit-learn 的普通最小二乘线性回归生成的𝛽₀和𝛽₁的值相比，Bambi 的贝叶斯线性回归生成的𝛽₀和𝛽₁的平均值如何？从寻找斑比创造的𝛽₀和𝛽₁的意义开始。

𝛽₀:的意义

0.43346411871369533

𝛽₁:的意义

7.403637994335104

它们非常接近 scikit-learn 生成的𝛽₀和𝛽₁的值！

绘制回归线

𝛽₀和𝛽₁.共有 4×1000 = 4000 个样本对于每一对𝛽₀和𝛽₁，我们可以得到𝑦.的一个估计

让我们用 10 对不同的𝛽₀和𝛽₁:来看 10 个预测

作者图片

酷！现在让我们使用𝛽₀ 的平均值和𝛽₁ 的平均值来得到𝑦(用黑线表示)的单一估计值。

作者图片

从上面的图中，我们可以看到黑线很好的呈现了所有的绿线。

所有预测的分布图

将单个预测分布与平均预测分布进行比较的另一种方法是将它们的分布相互重叠。

作者图片

从上面的图中我们可以看到，橙色线(平均预测分布)很好地代表了黑线(观察值)。

绘制预测的最高密度区间

对于𝑥 的特定值，最常见的预测是什么？让我们使用az.plot_hdi (HDI 代表最高密度区间)来解决这个问题。

作者图片

在上面的图中，

• 暗红色区间是给定𝑥值时 38%的预测值所在的区间
• 给定𝑥值，浅红色区间表示 68%的预测值

例如，如果数学成绩是 10，我们可以说葡萄牙语成绩的 38%的预测位于 10.5 到 12.8 之间。

通过创建贝叶斯线性回归模型，我们能够在给定𝑥.值的情况下创建大多数预测所在的区间很酷，不是吗？

结论

恭喜你！您刚刚学习了如何使用 Bambi 获得线性回归模型的系数分布，这使您能够获得可能预测的分布。

有一个分布将会给预测带来更多的洞察力。既然使用 Bambi 创建贝叶斯线性回归模型很容易，为什么不试一试呢？

本文的源代码可以在这里找到:

https://github.com/khuyentran1401/Data-science/blob/master/statistics/bayes_linear_regression/linear_regression.ipynb

我喜欢写一些基本的数据科学概念，并尝试不同的算法和数据科学工具。你可以通过 LinkedIn 和 Twitter 与我联系。

如果你想查看我写的所有文章的代码，请点击这里。在 Medium 上关注我，了解我的最新数据科学文章，例如:

参考

科尔森，W. (2018 年 4 月 20 日)。贝叶斯线性回归简介。中等。检索于 2021 年 12 月 29 日，来自https://towards data science . com/introduction-to-Bayesian-linear-regression-e66e 60791 ea 7

贝叶斯逻辑回归

原文：https://towardsdatascience.com/bayesian-logistic-regression-53df017ba90f?source=collection_archive---------9-----------------------

从零开始用朱莉娅语言

变化参数分布的模拟。图片作者。

如果您曾经搜索过评估模型准确性的评估指标，那么您很可能会发现许多不同的选项可供选择。从某种意义上说，准确性是预测的圣杯，所以机器学习社区花大量时间思考它一点也不奇怪。在一个越来越多的高风险决策被自动化的世界中，模型准确性实际上是一个非常有效的关注点。

但是，这种模型评估的方法看起来像是自动化决策制定的一种可靠而完整的方法吗？我们没有忘记什么吗？有些人会说，我们需要更多地关注模型的不确定性。无论你对你的模型进行了多少次交叉验证，优化的损失指标以及它的参数和预测仍然是固有的随机变量。仅仅关注预测的准确性而完全忽略不确定性会给自动决策系统带来错误的信心。因此，任何从数据中学习的可信的方法至少应该对其自身的不确定性保持透明。

我们如何估计模型参数和预测的不确定性？ Frequentist 不确定性量化方法通常涉及基于渐近假设或自举的封闭型解决方案(例如，参见此处的逻辑回归的)。在贝叶斯统计和机器学习中，我们关心的是对模型参数的后验分布建模。这种不确定性量化方法被称为贝叶斯推理，因为我们以贝叶斯方式处理模型参数:我们基于先验知识或信念对其分布做出假设，并根据新证据更新这些信念。频率主义者的方法避免了明确先前信念的需要，这在过去有时被认为是不科学的。然而，频率主义者的方法有自己的假设和缺陷(例如，参见 Murphy (2012))的讨论)。不深入讨论这个问题，现在让我们看看贝叶斯逻辑回归是如何自下而上实现的。

基本事实

在本帖中，我们将使用由二进制标签和相应的特征向量组成的合成玩具数据集。使用合成数据的好处是，我们可以控制生成数据的基本事实。特别地，我们将假设二元标签确实是由逻辑回归模型生成的。从混合高斯模型中生成特征。

为了给我们的例子增加一点生气，我们假设二进制标签根据猫和狗的身高和尾巴长度将样本分为猫和狗。下图显示了二维特征域中的合成数据。在下一节介绍贝叶斯逻辑回归之后，我们将使用这些数据来估计我们的模型。

地面真相标签。图片作者。

数学

在 Medium 上的文章中加入数学仍然有点麻烦，所以这里我们将完全依靠直觉，完全避免公式。使用 Julia 语言的好处之一是它允许使用 Unicode 字符，所以我们将在下面看到的代码看起来几乎和数学一样。如果你想看到完整的数学处理，以全面理解所有的细节，请随时查看我的网站上这篇文章的扩展版本。

问题设置

贝叶斯逻辑回归的起点是贝叶斯定理，它正式声明参数的后验分布与两个量的乘积成比例:给定参数的观察数据的可能性和参数的先验密度。应用到我们的上下文中，可以直观地理解如下:我们关于逻辑回归系数的后验信念由我们的先验信念和我们观察到的证据(即数据)形成。

在假设单个标签-特征对独立地分布和相同地分布的情况下，它们的联合似然性就是它们单个密度的乘积(伯努利)。关于参数的先验信念由我们决定。在实践中，它们可能来自以前的实验。这里我们将使用零均值球面高斯先验，原因将在下面进一步解释。

与线性回归不同，估计或最大化后验概率没有封闭形式的解析解，但幸运的是，精确的近似值确实存在(Murphy 2022)。一种最简单的方法叫做拉普拉斯近似法,实现起来非常简单，计算起来也非常高效。它依赖于以下观察结果:在高斯先验假设下，逻辑回归的后验概率也近似为高斯分布:特别是，该高斯分布以最大 T2 后验概率(MAP)估计为中心，其协方差矩阵等于在该模式下评估的逆 Hessian 矩阵。下面我们将看到如何估计地图和相应的黑森。

解决问题

实际上，我们不会直接最大化后验概率。相反，我们最小化负对数似然，这是等价的，并且更容易实现。下面的 Julia 代码显示了这个损失函数及其导数的实现。

免责声明 ❗️ 我应该提到这是我第一次用 Julia 编程，所以对于任何一个 Julia 的专业人士:请容忍我！非常乐意在评论中听到你的建议。

如您所见，负对数似然性等于两项之和:首先是(对数)伯努利分布之和(对应于数据的似然性)，其次是(对数)高斯分布(对应于我们先前的信念)。

拉普拉斯近似下贝叶斯逻辑回归的损失函数及其导数。

由于最小化这个损失函数是一个凸优化问题，我们有许多有效的算法可供选择来解决这个问题。有了 Hessian，使用二阶方法似乎是很自然的，因为包含关于损失函数曲率的信息通常会导致更快的收敛。这里我们将这样实现牛顿法:

Julia 语言中带 Arminjo 回溯的牛顿法。

后验预测

假设现在我们已经使用训练数据训练了贝叶斯逻辑回归模型作为我们的二元分类器，并且新的未标记样本到达。与任何二元分类器一样，我们可以通过简单地将新样本插入到我们的拟合模型中来预测丢失的标签。如果在训练阶段，我们发现模型达到了很好的精度，我们可以期待良好的样本外性能。但是因为我们仍然在处理一个随机变量的期望值，我们通常想知道这个预测有多嘈杂。

形式上，我们对后验预测分布感兴趣，它没有任何进一步的假设，是一个数学上难以处理的积分。它可以通过蒙特卡洛进行数值估计——通过简单地从后验分布中重复采样参数——或者通过使用所谓的**概率单位近似值。**后者利用了 sigmoid 函数可以由重新调整的标准高斯 cdf 很好地近似的发现(见下图)。以这种方式近似 sigmoid 函数允许我们导出后验预测的解析解。这种方法用于生成下一节中的结果。

概率单位近似的演示。图片作者。

估计值

下面的第一张图显示了在不同程度的先验不确定性下，高度和尾长系数的后验分布。红点表示无约束最大似然估计(MLE)。注意，随着先验不确定性趋于零，后验不确定性接近先验不确定性。这是很直观的，因为我们已经强加给自己，我们对先前的信念没有不确定性，因此再多的新证据也不能把我们推向任何方向。相反，对于非常大的先验不确定性，后验分布以无约束极大似然估计为中心:先验知识非常不确定，因此后验分布由数据的似然性决定。

不同先验不确定度σ下的后验分布。图片作者。

后验预测呢？情况类似:因为对于非常低水平的先验不确定性，后验完全由零均值先验支配，所以所有样本都被分类为 0.5(下图左上图)。随着我们逐渐增加先验的不确定性，预测后验越来越依赖于数据:预测标签的不确定性仅在没有样本填充的区域高。毫不奇怪，这种影响对于 MLE 来说是最强的，我们看到了一些过度拟合的证据。

不同先验不确定度σ下的预测后验分布。图片作者。

包扎

在这篇文章中，我们看到了贝叶斯逻辑回归是如何在 Julia 语言中实现的。模型参数的估计后验分布可用于量化系数和模型预测的不确定性。我认为对模型的不确定性保持透明是很重要的，以避免对估计值过于自信。

贝叶斯(概率)机器学习有更多的好处。了解我们的模型在输入域的什么地方表现出高度不确定性，例如，可以有助于标记数据:例如，参见 Gal、Islam 和 Ghahramani (2017)以及对图像数据的主动学习的有趣应用的后续工作。类似地，最近有一项工作在**算法追索权的背景下使用了后验预测的估计值(**舒特等人，2021)。关于算法追索权的简要介绍，请参见我的上一篇文章。

作为进一步阅读概率机器学习的重要参考，我强烈推荐 Murphy (2022)。这本书的电子版目前可以作为草稿免费获得。最后，如果你很好奇想详细看看完整的源代码，想亲自尝试一下代码，可以看看这个互动笔记本。

参考

克里斯托弗·毕晓普，2006 年。模式识别与机器学习。斯普林格。

加尔，亚林，里亚沙特伊斯兰和邹斌 Ghahramani。2017."图像数据的深度贝叶斯主动学习."在机器学习国际会议，1183–92。PMLR。

凯文·墨菲，2012 年。机器学习:概率视角。麻省理工出版社。

— -.2022.概率机器学习:简介。麻省理工出版社。

舒特，丽莎，奥斯卡·基，罗里·麦克·格拉斯，卢卡·科斯塔贝洛，波格丹一世·萨卡利安努，亚林·加尔，等。"通过隐含最小化认知和随机不确定性产生可解释的反事实解释."在国际人工智能和统计会议上，1756–64。PMLR。

全文发表于 2021 年 11 月 15 日 https://www.paltmeyer.comhttps://www.paltmeyer.com/post/bayesian-logistic-regression/。

通过 PyMC3 在 Python 中进行贝叶斯营销组合建模

原文：https://towardsdatascience.com/bayesian-marketing-mix-modeling-in-python-via-pymc3-7b2071f6001a?source=collection_archive---------0-----------------------

营销分析

一次性估计饱和度、残留和其他参数，包括它们的不确定性

由格雷格·拉科齐在 Unsplash 上拍摄的照片

在这篇文章中，我想结合我在早期文章中讨论的两个概念:贝叶斯建模和营销组合建模。因为你很有可能对这两个主题都不熟悉，所以让我给你一个快速的介绍和进一步的阅读。我会

1. 激发什么是营销组合建模，
2. 贝叶斯建模是什么，以及
3. 为什么将两者结合起来是有意义的。

然后，我将向您展示如何使用 PyMC3 进行实践。

准备工作

如果你是我文章的热心读者(谢谢！)，你可以跳过几节直接看代码。否则，请继续阅读。

营销组合建模

每个企业的一个基本问题是决定在哪些渠道上花费营销预算。你可以每天花 1000€在电视广告上，2000€在广播广告上，3000€在网页横幅广告上，跟着感觉走。但这有什么好处吗？

也许网页横幅频道已经在 T21 饱和了，在那里只花 1500€和 3000€一样好。这样你就可以节省 1500€，或者把它们投入到其他更好的销售渠道。

或者也许一些频道甚至有负的投资回报率——你在广告上花的每一欧元，你得到的回报不到一欧元。我们绝对不应该在这样一个频道上浪费太多的钱，至少从商业角度来看它没有战略重要性。

要回答这样的问题，你必须了解不同的媒体支出(电视、广播……)如何影响你的销售或其他感兴趣的 KPI。

在营销组合建模中，你从媒体支出数据集开始。它通常通过一些控制变量进行扩展，即关于任何可能影响目标 KPI 的信息，如假期、天气、足球锦标赛、停工、产品价格等等。不过，为了简洁起见，我们将省略控制变量。然后，当然，你需要一个你想要预测的 KPI。这通常是销售额、新客户数量等等。因此，典型的数据集可能如下所示:

图片由作者提供。

在我以前的文章中，我更详细地描述了动机以及如何进行营销组合建模。为了理解本文的其余部分，请在这里查看它们:

贝叶斯建模

许多估计器和模型源于一种最大似然方法。举个例子，假设你想估计一枚硬币正面朝上的概率 p 。你翻转 10 次，看到 8 个头，你得出什么结论？概率的自然估计是 p = 8 / 10 = 80%，这也是最大似然估计。你也可以计算一个置信区间，看看这个估计是否可靠，但是我们想走另一条路。

想象一下，我们想要加入一些关于概率 p 的先验知识。如果你从钱包里随机抽取了硬币，比如没有理由认为硬币有偏差，那么 p 应该不会离 50%太远，假设你不是魔术师。

使用贝叶斯建模，您可以将这种先验知识结合起来，最终得到对 p 的密度估计，即不是单个值，而是整个分布。这个分布可能在最大似然估计和先验之间的某处有一个峰值，可能是 65%。

红色曲线的模式是最大似然估计。图片由作者提供。

总而言之，贝叶斯建模是关于在先验知识和观察到的数据之间找到一个平衡点。在上图中，这意味着:没有任何数据，我们从蓝色曲线开始。这只是一种信念，一种直觉。然后，我们观察数据，告诉我们将蓝色曲线向红色曲线移动。我们以黄色混合曲线结束，它描绘了所谓的后验分布。

你可以在这里阅读更多关于动机的信息:

现在，理解理论是好的，但我们也必须能够应用它来完成事情。通常，我使用令人敬畏的 python 库 PyMC3 进行贝叶斯建模。你可以在这里看到它的作用:

为什么营销混合建模与贝叶斯？

您可以用许多超参数来定义营销组合模型:

• 饱和度
• 结转强度
• 结转长度
• …

然后，可以使用超参数优化方法来找到最优组合。这就是我在另一篇关于营销组合建模的文章中所做的， 用 Python 升级的营销组合建模。

这种方法很好，但是有一点我不喜欢:

超参数估计通常是不稳定的。

这意味着完全不同的超参数集可能产生同样好的模型。可能会有

• 型号 A 的电视残留强度为 0.4，电视饱和度为 0.8，以及
• 型号 B 电视结转强度 0.9，电视饱和度 0.5，

两者在测试装置上具有相同的 r 或 MAPE。从预测的角度来看，这两种模型是可以互换的，如果你保持在你目前所看到的营销支出范围内。

然而，使用模型 A 进行外推完全不同于使用模型 b 进行外推。这是一种相当不令人满意和麻烦的行为，因为外推是优化媒体预算时要做的事情。如果你过去一直每天花 0€到 1000€在电视广告上，为了优化，你必须知道当你花 5000€甚至 10000€时会发生什么。

你需要一个具有出色推断能力的模型。

通常，你不得不在两个以上的型号中选择。在这种情况下，您至少可以通过两种方式进行操作

• 你可以选择你创建的第一个模型，因为你甚至没有意识到这个问题。这种方法简单但危险。
• 您可以选择一个对您、一些领域专家或利益相关者来说正确的模型。对于一些人来说，这是可以的，但我不喜欢将产出预期放入模型中，因为下面的问题:

如果有人已经知道了答案，我为什么还要建立一个复制这个答案的模型呢？

从这里也可能有一些合理的方法，但是我现在想向你展示如何使用贝叶斯建模来规避这个问题。

建模开始！

首先，让我们获取数据集。

import pandas as pd

data = pd.read_csv(
  'https://raw.githubusercontent.com/Garve/datasets/4576d323bf2b66c906d5130d686245ad205505cf/mmm.csv',
  parse_dates=['Date'],
  index_col='Date'
)

X = data.drop(columns=['Sales'])
y = data['Sales']

然后，我们必须定义饱和度和结转功能，类似于上一篇文章。在 PyMC3 语言中，它可能是这样的:

import theano.tensor as tt

def saturate(x, a):
    return 1 - tt.exp(-a*x)

def carryover(x, strength, length=21):
    w = tt.as_tensor_variable(
        [tt.power(strength, i) for i in range(length)]
    )

    x_lags = tt.stack(
        [tt.concatenate([
            tt.zeros(i),
            x[:x.shape[0]-i]
        ]) for i in range(length)]
    )

    return tt.dot(w, x_lags)

饱和度函数应该很容易掌握。然而，结转有点麻烦。基本上，您可以将遗留转换表示为矩阵向量乘法。你只需要先组装矩阵x_lags和向量w就可以了。作为一个例子，我们可以转换输入向量 x = ( x ₁， x ₂， x ₃， x ₄)，结转长度为 3 via

图片由作者提供。

上面代码中的carryover函数就是这样做的。有了这些功能，我们终于可以开始建模了。

import pymc3 as pm

with pm.Model() as mmm:
    channel_contributions = []

    for channel in X.columns:
        coef = pm.Exponential(f'coef_{channel}', lam=0.0001)
        sat = pm.Exponential(f'sat_{channel}', lam=1)
        car = pm.Beta(f'car_{channel}', alpha=2, beta=2)

        channel_data = X[channel].values
        channel_contribution = pm.Deterministic(
            f'contribution_{channel}',
            coef * saturate(
                carryover(
                    channel_data,
                    car
                ),
                sat
            )
        )

        channel_contributions.append(channel_contribution)

    base = pm.Exponential('base', lam=0.0001)
    noise = pm.Exponential('noise', lam=0.0001)

    sales = pm.Normal(
        'sales',
        mu=sum(channel_contributions) + base,
        sigma=noise,
        observed=y
    )

    trace = pm.sample(return_inferencedata=True, tune=3000)

我们可以看到所有参数(不再有超参数！)用粗体标出。它是回归系数、饱和强度、残留强度、基线和噪音。

请注意，我没有考虑结转长度，而是将其设置为 21。这是因为我还不知道如何在 PyMC3 中用可变维度在carrover函数中创建矩阵和向量。但通常，我会用泊松随机变量来表示遗留长度。如果你知道如何正确地做这件事，请给我留言！😉也许现在也是给 pyro 一个机会的好时机，PyMC3 这样的另一种概率编程语言。

模型输出的分析

之后，我们可以看看通常的图片。让我们从后验分布开始。执行

import arviz as az

az.plot_posterior(
    trace,
    var_names=['~contribution'],
    filter_vars='like'
)

生产

图片由作者提供。

在这里，你可以看到所有参数的后验概率。它们都有很好的单峰形状。您还可以通过以下方式探索变量对的行为

az.plot_joint(
    trace,
    var_names=['coef_TV', 'sat_TV'],
)

图片由作者提供。

这里你可以看到饱和强度和回归系数不是独立的，而是负相关的:系数越高，饱和参数越低。这是有意义的，因为较高的系数可以补偿较慢增加的饱和曲线(=较低的sat_TV),反之亦然。

让我们看另一个例子:

az.plot_joint(
    trace,
    var_names=['car_TV', 'sat_TV'],
)

图片由作者提供。

在这里，我们可以看到为什么超参数优化可能会有问题。这张图片中的每个点都是一个潜在的模型，你可以通过超参数优化找到它。

对于一个真正唯一的最佳模型，我们更愿意看到一个点云紧紧的集中在一个点上(car_TV_true，sat_TV_true)。然而，在这里，我们看到 TV 残留强度可以具有在 0.4 和 0.5 之间的合理值，这取决于饱和度参数。

在我们下结论之前，我们还可以检查一下这个模型是否好。执行

import matplotlib.pyplot as plt

with mmm:
    posterior = pm.sample_posterior_predictive(trace)

means = posterior['sales'].mean(0)
stds = posterior['sales'].std(0)

plt.figure(figsize=(20, 8))
plt.plot(y.values, linewidth=2, c='r', label='Observations')
plt.plot(means, linewidth=1, c='b', label='Mean prediction')
plt.fill_between(np.arange(len(y)), means - 2*stds, means + 2*stds, alpha=0.33)
plt.legend()

给了我们

图片由作者提供。

所以看起来这个模型找到了一些有用的东西。在这里，我不会进一步详细讨论如何评估模型的性能，我们可以在将来这样做。

渠道贡献

到目前为止，我们已经处理了分布，但是对于我们最喜欢的渠道贡献图，让我们再次以单个值结束。由于我们在 PyMC3 代码中引入了一些通道贡献变量，现在我们可以使用一个简短的compute_mean函数轻松提取它们。

def compute_mean(trace, channel):
    return (trace
            .posterior[f'contribution_{channel}']
            .values
            .reshape(4000, 200)
            .mean(0)
           )

channels = ['Banners', 'Radio', 'TV']
unadj_contributions = pd.DataFrame(
    {'Base': trace.posterior['base'].values.mean()},
    index=X.index
)

for channel in channels:
    unadj_contributions[channel] = compute_mean(trace, channel)

adj_contributions = (unadj_contributions
                     .div(unadj_contributions.sum(axis=1), axis=0)
                     .mul(y, axis=0)
                    )

ax = (adj_contributions
      .plot.area(
          figsize=(16, 10),
          linewidth=1,
          title='Predicted Sales and Breakdown',
          ylabel='Sales',
          xlabel='Date'
      )
     )

handles, labels = ax.get_legend_handles_labels()
ax.legend(
    handles[::-1], labels[::-1],
    title='Channels', loc="center left",
    bbox_to_anchor=(1.01, 0.5)
)

图片由作者提供。

看起来不错！

总结与展望

在本文中，我们讨论了使用最大似然法的营销组合建模可能会因为超参数估计而变得困难。可能有许多模型在测试集上表现良好，但具有完全不同的外推行为——这些(超)参数非常不稳定。

然而，适当的外推是优化的关键。因此，我们设计了一个基本的营销组合模型，在贝叶斯设置中具有饱和和结转效应。这被证明是有用的，因为它一次估计所有的参数，给我们更稳定的参数估计。

然后，我们使用 PyMC3 实现了它，并再次创建了一个漂亮的贡献图表，以查看渠道支出对销售额的贡献。

模型扩展

我们现在可以采用我们已经创建的贝叶斯模型，并进一步扩展它。例如，我们可以引入随时间变化的参数。如果两年前电视遗留强度开始为 0.8，但随着时间的推移慢慢下降到 0.5，这可能是有用的，这是概念漂移的一个例子。我们可以使用高斯随机游走，轻松地对此建模，如我关于滚动回归的另一篇文章所示:

同样，我们可以对变化的基线进行建模。到目前为止，我们已经把基线当作一个单一的数字，它在整个训练期间是固定的。但也许随着时间的推移，由于我们的战略广告努力，它也增加了。这是使用我们旧的超参数优化最大似然模型很难解决的问题。

我希望你今天学到了新的、有趣的、有用的东西。感谢阅读！

作为最后一点，如果你

1. 想支持我多写点机器学习和
2. 无论如何都要计划获得中等订阅量，

为什么不做 通过这个链接 ？这将对我帮助很大！😊

说实话，给你的价格没有变化，但大约一半的订阅费直接归我。

非常感谢，如果你考虑支持我的话！

有问题就在 LinkedIn 上写我！

网球运动员排名的贝叶斯模型

原文：https://towardsdatascience.com/bayesian-modelling-for-tennis-player-ranking-d151aef789a7?source=collection_archive---------27-----------------------

贝叶斯模型如何帮助我们对网球比赛结果建模并创建运动员排名？

图片来自 Pixabay 上的 Bessi

我的前两篇文章是关于马尔可夫链蒙特卡罗和变分推理背后的理论，它为贝叶斯建模奠定了基础。在本文中，我们将使用贝叶斯模型和 MCMC 来预测网球比赛的结果，并在这样做的过程中，创建一个网球运动员的排名。

数据如何？

本项目使用的所有数据均来自此处的。我选择使用 2019 年的数据，因为 2020 年对体育来说是一个怪异的年份(&其他一切！)

2019 年的数据包含 2611 场男子单打比赛，来自所有 ATP(网球职业选手协会)锦标赛。对于这个项目，我们将只使用数据的Winner、Loser和Surface列。

本项目所需的数据处理如下:

1. 将玩家姓名映射到玩家 ID
2. 将曲面映射到曲面 ID
3. 将Winner和Loser列改为player1和player2
4. 随机洗牌，使player1不总是赢家(我们不希望任何模型知道player1有优势)。
5. 创建一个名为result_category的新列，指示player1是否赢得了游戏。

这给我们留下了一个看起来像这样的数据帧:

经过处理的 ATP 网球数据，图片由作者提供

建模理论

这个问题的建模方法围绕着一个想法，即一个网球运动员有一个潜在的技能分数，决定了他们赢得比赛的能力。如果一个玩家的分数比另一个玩家高，他们应该更有可能赢。此外，一个玩家应该有 3 个不同的奖励，取决于他们玩的表面(草地，粘土或硬)。这应该模拟玩家在不同表面上的真实可变性。

值得注意的是，由于我们将此视为贝叶斯建模问题，一切都变成了概率，所有技能和表面奖励都带有不确定性，这应该模拟不同玩家的真实方差。

实践建模

有很多不同的包可以用来解决这样的问题，包括PyMC3和pyro.，我选择使用PyStan，它是一个包装Stan编程语言的包装器。

我不会深入讨论如何使用 Stan 的细节，因为 Stan 手册比我能做的要好得多，但是我会试着解释一些基本的东西。有 3 个主要模块定义一个模型，data、parameters和model(可选地，你也可以包括generated quantities)。

1. data定义所有的模型特征和目标(如果有的话)
2. parameters定义所有你想要取样的变量
3. model定义参数先验和模型似然
4. generated quantities在采样后运行，用于执行模型推断

主要关注的是model块，我们在这里定义先验和似然性。

似然性表明给定的网球比赛分布为二元伯努利随机变量，其中参数 p 由运动员能力差异的 sigmoid 确定。

来自标准文档的伯努利逻辑函数。注意 logit^-1 和乙状结肠一样

玩家能力是两件事的函数，基础技能分数player_ability和表面奖励surface_bonus。注意选择不使用超先验参数化player_ability和surface_bonus的平均值。这是因为可能性仅取决于player_ability和surface_bonus中的差异，因此任何关于平均值的先验只会导致约束较差的问题。

请注意，我们确实允许每个表面有单独的标准偏差，即一些表面会导致或多或少的性能差异。最后值得注意的建模细节是两个标准差的先验选择。两者都被选作uniform(0,10)，这是一个完全未知的先验，目的是不在模型上强加太多的先验信息。我们同样可以选择像逆伽马先验这样的东西，但是我在采样时似乎收敛得稍微差一些。

结果

因此，如果我们对大约 90%的游戏(随机选择)进行训练，并对最后 10%进行验证，我们在测试集上获得的二进制准确率大约为 70%，在维持集上获得的二进制准确率为 65%。这并不坏，因为我们给模型的唯一信息是谁玩了，结果是什么，他们在什么表面上玩。

但是请记住，贝叶斯模型允许对我们预测的不确定性进行测量。如果我们根据预测样本的标准偏差来设定预测的阈值，我们可以建立这样的曲线…

σ与精度图，图片由作者提供

因此，通过设定大约 0.49 的标准偏差阈值，我们可以将准确率提高到近 90%！

玩家排名

我们的模型要求我们学习每个玩家的分数，以便预测谁将赢得给定的网球比赛。所以我们的模式还免费给了我们一个玩家排名！

模特选手排名，作者图片

如果把排名和 2019 年底开始的 ATP 排名对比一下…

2019 年底 ATP 排名

我们以正确的顺序获得了前 6 名，并且我们正确地找到了 ATP 前 10 名中的 8/10！对于几个小时的造型来说还不错！

我们学到了什么？

在本文中，我们在 Stan 中构建了一个贝叶斯模型来预测网球比赛的结果。我们还研究了如何使用该模型为球员创建一个与 ATP 排名相当一致的球员排名！使用 Stan 既直观又有趣。我很期待看到我们还能用它做些什么！

贝叶斯网络:分析酒店客户数据

原文：https://towardsdatascience.com/bayesian-networks-analysing-hotel-customer-data-cd67dfd6fd43?source=collection_archive---------16-----------------------

用 bnlearn 实现概率建模

来源:图片由 geralt 从 Pixabay 拍摄

在检查不同变量之间的依赖关系时，贝叶斯网络是一个非常直观的工具。

具体来说， DAG (或有向无环图)是允许我们表示给定变量集之间的条件概率。

使用 Python 中的 bnlearn 库，让我们看看这样的分析如何让我们收集关于酒店客户不同属性之间的依赖性的见解。

背景

本分析中使用的酒店预订数据来自 Antonio、Almeida 和 Nunes (2019)。

采用不同酒店客户属性的数据集，如客户是否取消、他们的原籍国、市场细分等，从数据集中随机抽取 1000 个观察值。

来源:Jupyter 笔记本输出，Antonio、Almeida 和 Nunes 的数据集(2019 年)

数据处理和 DAG 生成

我们将首先把数据转换成一个热数据矩阵。该函数通过自动检测相关数据类型并在必要时实现一键编码来工作:

dfhot, dfnum = bn.df2onehot(df)

现在，我们将生成 DAG。

首先，将根节点设置为is cancelled，即客户是否取消了酒店预订。在这种情况下，我们将它设置为根节点——因为我们对研究不同的客户属性如何影响取消行为特别感兴趣。

此外，我们还选择放弃某些变量，如提前期，因为我们主要对分析整个数据集的分类特征感兴趣，并根据每个特征中包含的不同类别分配取消概率。

与数据集中其他要素紧密关联的要素将被丢弃。例如， ReservedRoomType 与 AssignedRoomType 紧密相关，因此决定放弃前者。

DAG = bn.structure_learning.fit(dfnum, methodtype='cl', black_list=['LeadTime', 'Country', 'DistributionChannel','ReservedRoomType'], root_node='IsCanceled', bw_list_method='nodes')G = bn.plot(DAG)

以下是生成的 DAG:

来源:Jupyter 笔记本输出

如上所述，我们现在可以看到变量之间依赖关系的网状图。

例如，取消似乎证明了一个很强的条件概率，其中市场细分和需要停车位。换句话说，某个细分市场的客户比其他人更有可能取消预订，停车位的可用性似乎也对客户是否选择取消预订有很大影响。

反过来， MarketSegment 变量似乎显示出与 AssignedRoomType 和 DepositType 等的显著相关性。例如，某些细分市场可能比其他市场更愿意提前支付押金，或者某些细分市场可能为某个房间类别支付更多费用。

当然，这是表面层次的分析，因为虽然 DAG 的可视化可以向我们显示不同变量之间存在依赖关系，但它们无法解释这些依赖关系的性质或它们存在的原因。

从这个角度来看，在进行进一步分析之前，对不同变量之间的关系有一个总体的了解是非常有用的。

进一步分析

在更详细地生成 DAG 的图之后，我们现在可以更详细地分析参数。

简单看一下另一个例子

让我们考虑一个不同的例子 bnlearn 库提供的洒水喷头例子。

来源:Jupyter 笔记本输出

当分析这里的相关性时，我们可以看到雨和云之间存在关系，无论洒水器是开着还是关着，以及湿草的发生率。

可以在变量之间计算 CPD(或条件概率分布)。

例如，这里是生成的 Rain 变量的 CPD。

[bnlearn] >CPD of Rain:
+---------+--------------------+---------------------+
| Cloudy  | Cloudy(0)          | Cloudy(1)           |
+---------+--------------------+---------------------+
| Rain(0) | 0.6518218623481782 | 0.33695652173913043 |
+---------+--------------------+---------------------+
| Rain(1) | 0.3481781376518219 | 0.6630434782608695  |
+---------+--------------------+---------------------+

如果是阴天，看不到下雨的概率超过 65%。在这种情况下，观测到降雨的概率不到 35%。

另一方面，如果有云，那么观察到下雨的概率超过 66%，观察到不下雨的概率略高于 33%。

回到我们酒店的例子

既然我们已经考虑了一个更简单的例子，那么让我们来看看在酒店数据中生成的一些 CPD。

>>> model = bn.parameter_learning.fit(DAG, dfnum)

看到取消本身的概率刚刚超过 38%。

[bnlearn] >CPD of IsCanceled:
+---------------+-------+
| IsCanceled(0) | 0.611 |
+---------------+-------+
| IsCanceled(1) | 0.389 |
+---------------+-------+

现在，让我们考虑不同变量对取消概率的影响。

所需停车位

来源:Jupyter 笔记本输出

我们可以看到，如果没有所需的停车位可用，那么客户取消的可能性会大大增加。

细分市场

来源:Jupyter 笔记本输出

从该 CPD 中，我们还可以看到某些细分市场显示出比其他市场更高的取消概率。

现在，让我们假设我们希望考虑选择变量组合对取消概率的影响。

我们以餐食变量为例。如果酒店不满足顾客的饮食要求，任何一位顾客取消预订的可能性有多大？

>>> query = bn.inference.fit(model, variables=['IsCanceled'], evidence={'Meal':True})
>>> print(query)
>>> print(query.df)+----+--------------+---------+
|    |   IsCanceled |       p |
+====+==============+=========+
|  0 |            0 | 0.62191 |
+----+--------------+---------+
|  1 |            1 | 0.37809 |
+----+--------------+---------+

那么，如果我们也考虑所需的停车位呢？

>>> query = bn.inference.fit(model, variables=['IsCanceled'], evidence={'RequiredCarParkingSpaces':True, 'Meal':True})
>>> print(query)
>>> print(query.df)+----+--------------+---------+
|    |   IsCanceled |       p |
+====+==============+=========+
|  0 |            0 | 0.66197 |
+----+--------------+---------+
|  1 |            1 | 0.33803 |
+----+--------------+---------+

我们可以看到取消的概率下降到 33%多一点。当去除用餐变量时，概率刚刚超过 34%。从这个角度来看，分析表明，一家酒店如果有足够的停车位供有需求的顾客使用，那么取消预订的可能性就会大大降低。

通过这种方式， bnlearn 库在提供变量之间存在的条件概率的高层次概述以及对那些依赖关系的本质的更详细分析方面非常有用。

结论

在本文中，您已经看到:

• 贝叶斯网络如何运作
• 如何绘制和解释定向丙烯图(DAG)
• 用 bnlearn 计算条件概率

非常感谢您的宝贵时间，非常感谢您的任何问题或反馈。你可以在这里找到原始文章以及该项目的 GitHub 文件夹的链接。

参考

• 安东尼奥、阿尔梅达、努内斯(2019)。酒店预订需求数据集。
• erdogant.github.io: bnlearn 文档。

免责声明:本文是在“原样”的基础上编写的，没有任何担保。它旨在提供数据科学概念的概述，不应被解释为专业建议。本文中的发现和解释是作者的发现和解释，不被本文中提到的任何第三方认可或隶属于任何第三方。作者与本文提及的任何第三方无任何关系。

贝叶斯神经网络

原文：https://towardsdatascience.com/bayesian-neural-network-7041dd09f2cc?source=collection_archive---------13-----------------------

图 1:示意图显示了这项工作中使用的贝叶斯神经网络的架构。该网络有一个包含八个参数的输入层、一个包含十二个节点的隐藏层和一个包含一个节点的输出层。输入层和隐藏层之间的权重由正态分布的 w0_ij 定义。I，j 是节点输入和隐藏层节点索引。类似地，w1_jk 是隐藏层和输出层之间权重的正态分布。在输出节点，网络产生介于 0 和 1 之间的预测分数分布。图摘自我在 arxiv 上发表的一篇文章:https://arxiv.org/abs/1911.09660

在传统的神经网络中，权重被指定为单个值或点估计值，而在 BNN 中，权重被视为概率分布。网络权重的这些概率分布用于估计权重和预测的不确定性。图 1 显示了重量呈正态分布的 BNN 的示意图。使用贝叶斯定理计算权重的后验值，如下所示:

其中 X 是数据，P(X|W)是观察到 X 的可能性，给定权重(W)，P(W)是权重的先验信念，分母 P(X)是数据的概率，也称为证据。它需要对所有可能的权重值进行积分，如下所示:

综合证据中所有的不确定权重，很难找到一个封闭形式的解析解。因此，基于模拟或数值的替代方法，如蒙特卡罗马尔可夫链(MCMC)，变分推理(VI)被考虑。MCMC 抽样已经成为现代贝叶斯统计中重要的推理方法。科学家们在许多应用中广泛研究和应用。然而，该技术对于大型数据集或复杂模型来说很慢。另一方面，变分推断(VI)比其他方法更快。它还被应用于解决许多大规模计算昂贵的神经科学和计算机视觉问题。

在 VI 中，考虑一个新的分布 Q(W|θ),它近似于真实的后验 P(W|X)。Q(W|θ)通过 W 进行参数化，VI 通过优化找到使两个分布的差异最小化的正确集合:

在上面的等式-3 中，KL 或 kull back–lei bler 散度是真实分布和近似分布之间相似性(相对熵)的非对称和信息论度量。Q(W|θ)和 P(W|X)之间的 KL 散度定义为:

用等式-1 代替 P(W|X ),我们得到:

取 Q(W|θ)的期望值，我们得到:

上面的等式仍然显示了 logP(X)的依赖性，使得很难被 KL 计算。因此，通过添加具有负 KL 散度的 logP(X) 来导出备选目标函数。 logP(X) 是关于 Q(W|θ)的常数。新函数称为下界证据(ELBO ),表示为:

第一项称为似然性，第二项是变分分布和先验权重分布之间的负 KL 散度。因此，ELBO 在可能性和先验之间进行平衡。可以使用不同的优化算法(如梯度下降)来优化 ELBO 目标函数以最小化 KL 散度。

想要在你的领域找到对机器学习和图论感兴趣的人吗？

我将通过分享关于 xoolooloo 的有趣信息来结束这篇文章。它是一个基于位置的搜索引擎，使用相似和多重兴趣来查找当地人。例如，如果你读了这篇文章，你肯定对图论、机器学习感兴趣。因此，你可以在你的地区找到有这些兴趣的人；去看看 www.xoolooloo.com 的

非常感谢您的阅读。我希望你喜欢这篇文章。完整代码和相关数据可在 Github 上找到。文章链接:https://arxiv.org/abs/1911.09660。我很想收到你的来信。你可以联系我:

Email: sabbers@gmail.com
LinkedIn: [https://www.linkedin.com/in/sabber-ahamed/](https://www.linkedin.com/in/sabber-ahamed/)
Github: [https://github.com/msahamed](https://github.com/msahamed)
Medium: [https://medium.com/@sabber/](https://medium.com/@sabber/)

贝叶斯优化:一步一步的方法

原文：https://towardsdatascience.com/bayesian-optimization-a-step-by-step-approach-a1cb678dd2ec?source=collection_archive---------3-----------------------

统计、概率、优化

用统计细节解释贝叶斯优化

M. B. M. 在 Unsplash 上拍摄的照片

在许多现实生活的分析用例中，优化函数是非常重要的。所谓优化，我们的意思是，用某一组参数组合找到目标函数的最大值或最小值。找出最小值或最大值以及参数应该是目标。在本文中，我们将讨论用贝叶斯方法优化一个未知的代价函数的基础。

昂贵的黑盒函数优化——微积分

从基础微积分我们知道，要找到一个函数的最大值或最小值，我们需要求解 x 的导数方程，如下所示:

方程的根将是参数 x 的最优值，这又给出了函数 f(x)的最优值。只要知道函数 f(x)的完整代数形式，就很简单。但是，在许多现实生活场景中，事情并非如此简单。如果它是一个黑盒，那么你将只知道任何 f(x)的输出&输入值，但是你将不知道完整的形式。所以，分析上你找不到导数。除此之外，调用该函数的成本可能非常高。考虑如下两个用例:

1.找出神经网络的最佳超参数组合

考虑一个解决分类问题的大型复杂神经网络。超参数可以是网络中隐藏单元的数量、隐藏层的数量等。这些之间的关系可以被认为是一个假设的函数，该函数将超参数作为输入，并将分类精度作为输出。当然，你不知道这个函数的实际代数形式。找到最佳组合的一种方法是通过反复训练网络来尝试各种随机组合。问题的根源就在于此。我们负担不起重复培训，因为这是一个非常庞大和复杂的网络，培训需要大量的时间和资源。贝叶斯优化可以在这里有所帮助。

2.考古遗址的挖掘——寻找最佳“挖掘点”

不仅对于软件(如神经网络案例)，贝叶斯优化也有助于克服物理世界的挑战。在一个考古现场，专家们想到的主要问题是:“在哪里挖掘？”。不用说，从人力、金钱和时间的角度来看，挖掘是一个昂贵的“黑箱”操作。这个函数可以被认为是返回了一个站点中可用资源的列表，参数可以是位置细节，以及其他一些特定领域的东西。找到那个位置是一个挑战，贝叶斯方法可以解决这个问题。

在所有情况下，都需要进行一些函数的初始运行来进行估计。考虑以下功能:

它是“白盒”(因为导数可以通过分析找到)，乍看起来成本并不高，但是为了方便起见，假设它是成本很高的&本质上是“黑盒”，需要大量时间来返回输出。让我们进行一些模拟，

x = np.random.randn(5,2)
y = costly_function(x)
pd.DataFrame(data={'y':y, 'x0':x[:,0], 'x1':x[:,1]})

结果呢，

该函数有两个参数 x0 和 x1。对于上面这种代价很高的函数，我们只能调用几次。可能发生的情况是，初始的函数结果&输入可能在一个文件或数据库中给我们。然后，我们必须找出最佳(最小或最大)值和与之相关的参数组合(本例中为 x0 & x1)。我们永远不会知道实际的代数形式，导数的解析形式，因此必须进行数值优化。

在贝叶斯优化中，输入/输出组合的初始集合通常如上所述给出，或者可以从函数中生成。对于上面讨论的两个用例，可以如下实现:

1. 神经网络在不同的超参数组合上被训练多次，并且精确度被捕获和存储。该集合可以用作初始数据点。
2. 在考古挖掘作业中，可以进行几次初步挖掘来收集有关遗址的信息。它作为初始数据点。

简而言之，这是一个约束优化，它解决了下面给出两个问题:

I)找出以数值方式给出黑盒函数最优值的最优参数，因为无法找到解析导数。

ii)保持整个过程中函数调用的数量尽可能的少，因为这是非常昂贵的。(除了最初的几次运行)

贝叶斯优化术语

贝叶斯方法基于“黑盒”函数的统计建模和参数空间的智能探索。了解几个术语是很重要的。

1.代理模型

它是“黑箱”功能的统计/概率模型。它是后者的代理。为了试验不同的参数，该模型用于模拟函数输出，而不是调用实际的高成本函数。高斯过程回归(它是多变量高斯随机过程)被用作贝叶斯优化中的“代理”。

2.采集功能

它是一个度量函数，决定哪个参数值可以从函数中返回最佳值。它有许多变体。我们将与“预期改进”一起工作。

3.勘探 vs 开采

典型的策略是在参数空间中的局部和全局最优值之间进行补偿。考虑下面给出的图表:

https://www.shutterstock.com/image-vector/local-global-maxima-minima-1961848822

在进行参数空间探索时，可以在函数具有高值或低值的地方找到许多这样的局部最优数据点。但是，这一过程不应就此停止，因为在其他一些领域可能有更多的最佳值。它被称为“探险”。另一方面，也应该重视从函数中持续返回最优(高或低)值的点。就是“剥削”。所以，两者都有一定的意义。这是一个微不足道的决定，“什么时候探索不同位置的更优数据点，或者什么时候利用同一个方向的&”。这是贝叶斯优化在参数空间中击败传统随机搜索或网格搜索方法的领域，因为它采取了中间立场。它有助于通过少量的实际函数调用更快地实现目标。另外两种方法是完全忽略这一事实的盲目搜索。搜索必须非常精确&“切中要害”以降低成本。贝叶斯方法很好地解决了这个问题。

简而言之，采集函数使用“探索对开发”策略，以迭代方式决定最佳参数搜索。在这些迭代中，代理模型有助于获得函数的模拟输出。任何贝叶斯方法都是基于“先验/后验”的概念。前面提到的函数的初始运行被用作起始点或“先验”，并且在每次迭代中，这些“先验”被“后验”数据点丰富。几次迭代后，达到最佳数据点，整个过程在此停止。接下来，我们将看到所有这些技术的实际应用。

逐步实施

我们将使用前面部分声明的带有两个参数的同一个“成本函数”作为要优化的目标函数，即，我们需要在当前情况下最大化它。假设该函数是“黑盒”,很少运行的输出会给我们。我们将定义一个类“BayesianOptimizer ”,并逐步声明它的功能和属性。

让我们首先创建类的结构并声明它的属性。

“gauss_pr”是前面提到的代理模型。“x_init”和“y_init”通常是初始运行集合中函数的参数和输出。您将及时看到其他属性的用法。

现在，我们将在同一个类中创建一个名为“_get_expected_improvement”的函数，这是贝叶斯方法的核心。就是前面说的想要的采集函数。在开始之前，让我们先讨论一些理论。

每当尝试一个新的数据点时，我们都需要计算一个称为“预期改进”或“EI”的指标来给出该数据点的权重。公式由下式给出:

在哪里

这乍一看有点复杂。“mu(f(x))”是根据高斯过程对新数据点 x 的预测(即“y”值)。“max{f(x)}”是根据当前阶段的整个先验列表的预测的最大值(也是从高斯过程获得的)。“sigma(f(x))”通常是新数据点 x 的预测标准偏差。“mu(f(x))”和“max{f(x)}”之间的差异只是为了检查搜索过程的改进。较高的值表示新数据点正在从函数返回一个高值，该值比迄今为止获得的最大值“显著”高。这个“显著性是通过将差值乘以累积概率密度得到的。因此，它给出了一个预期的或整体的"表示改进，当然这是"开发部分。这个值还被一个神奇的因子“sigma(f(x))”放大了。它给出了“EI”度量的不确定性，并且是处理“探索部分的秘密。“mu(f(x))”和“sigma(f(x))”很好地平衡了彼此的影响，因此 EI(x)采取了中间立场。

现在，我们将看到这个函数的实现，

它避免了实际的函数调用，并使用高斯过程作为代理。因此，具有不同点的试验通过代理高斯过程发生，而不是实际函数。这就是“降低成本”的秘诀之一。

这个高斯过程是如何迭代建立的，我们后面会看到。

接下来，这个获取函数用于计算随机选择的数据点的邻域中的“EI”度量，在数值上，它通过数值导数计算而被最大化。不要在这里混淆！！我们刚才又说了函数的“数值导数”。但是，它不是目标成本函数。它是关于最大化采集函数。只要想一想！！我们将只对给出“EI”最大值的数据点感兴趣，即从当前最大值给出目标函数的最大改进(“y”值)。

当然，这将为我们继续搜索参数空间提供正确的方向，避免不必要的盲目探索。这是“降低成本”的第二个秘诀。

接下来我们将看到实现，

尝试了一批(由“batch_size”提到的)随机生成的数据点，并且对于每一个数据点，采集函数被最大化(实际上，负函数被最小化以匹配“scipy”库支持)。将函数的负值最小化等价于将其最大化)。同样，从迭代中取“EI”的最大值，并返回相应的参数 x 值。从一个起点开始数值最大化，再从结果中取最大值，实际上给了参数 x 一个正确的方向，有适量的平均值和不确定性。它是广度搜索(不确定性)和深度搜索(均值)之间的折衷，从两者中获得最佳效果。

现在，下一部分是实际工作。您看到有一个高斯过程回归器作为代理模型。它用于在“先前”数据点的基础上迭代构建模型。“优化”功能可以完成所有这些工作。如下图所示。

上述步骤可以总结如下:

I)从初始“先前”数据点中找出最大值和相应的参数 x

ii)用初始“先前”数据点构建高斯过程。高斯过程使用最大似然估计来寻找参数之间的正确联合分布。

iii)使用采集函数获得下一个最佳参数 x。如前所述，在这个步骤中，在高斯模型的帮助下对不同的数据点进行试错，而不需要调用实际的“昂贵的”目标函数。

iv)使用真实目标函数获得参数值 x 的“y”值。它是“后验”数据点。

v)用“后验”数据更新“先验”,根据“先验”更新当前最大值，然后再次转到步骤 ii)

vi)最后，返回当前最大“y”值和相应的参数 x。

您会注意到，我们还进行了一些距离计算，并捕获了当前的最佳样本(除了当前的最大值 y 和相应的 x 之外，什么也没有)。从函数“_get_next_probable_point”计算两个连续的可能的下一个 x 值之间的距离。它显示了算法的进度。你会在结果中看到这一点。

结果

我们将使用两个二维样本参数作为初始数据点。我们可以使用 n 维，因为“costing _ function”足够通用，可以处理这种情况。

sample_x = np.array([[8,1],[6.2,5.3]])
sample_y = costly_function(sample_x)

现在，是时候使用“BayesianOptimizer”类了，

bopt = BayesianOptimizer(target_func=costly_function, x_init=sample_x, y_init=sample_y, n_iter=200, scale=10, batch_size=30)
bopt.optimize()

它触发邻域大小为 30 的搜索，并进行 200 次迭代。结果如下所示:

因此，对于参数 x0=1.92 和 x1=3.62，代价函数的最大值是 3.7。你可能已经注意到，在代价函数的定义中，我们引入了一个随机噪声，只是为了让它的导数难以计算。

记住，正如前面多次说过的，你永远也不会得到“高成本函数”的主体或定义。这里，为了便于理解，我们只创建了一个虚拟函数。但是，实际上，情况并非如此。样本 x & y 可能作为数据集提供给你，或者会有一些公共/私有托管的“黑盒”API，可以调用它来获得任何 x 的 y 值。

现在，我们将从计算中得到一些其他的结果，

可以画出距离，

pd.DataFrame(bopt.distances_).plot()

注意，两个连续的可能的下一个 x 点之间的距离随着迭代而减小。这是意料之中的。高斯模型在每次迭代后变得更加成熟，其预测变得更加完美，从而产生精确的“EI”值。最终，在每次迭代之后，它更接近最优 x，因此距离开始减小。

现在，我们将看到“y”值在迭代中是如何变化的，

bopt.best_samples_['y'].plot()

注意，y 是逐步达到最大值的。每次迭代都产生对最大值“y”的更好估计。

现在，同样的图计算“EI”

bopt.best_samples_['ei'].plot()

“EI”正在如预期的那样下降。想想吧！！。每一次迭代之后，都会产生一个更好的 y，所以这意味着当我们向最优的方向前进时，获得更大改进的机会就更少了。

有趣的事实

一个大问题“你如何决定 n_iter？”改变它的值肯定会影响我们对最大值的估计。这里需要注意的一点是，

n_iter 等于在除初始化部分之外的整个优化过程中调用实际代价函数的次数

所以，当然它的价值取决于我们准备承担多少成本。

考虑神经网络模型的用例 1。如果有一个云环境，它对深度神经网络的每次“训练”运行向其客户收费，并且存在预算约束，则必须相应地设置“n_iter ”(因为在这种情况下，目标函数是网络的训练&获得准确性)。

贝叶斯方法试图通过减少实际调用来给出函数的估计，因此在某些情况下其准确性可能不如 RandomSearch 或 GridSearch。

当成本比非常小的精确度更重要时，贝叶斯优化是有用的。

源代码可以在这里找到，

https://github.com/avisheknag17/public_ml_models/blob/master/bayesian_optimization/bayesian_optimization.ipynb