如何选择咱们的损失函数 - 链接统计推断和问题所在领域的桥梁

1. 从几率预测到决策函数之间的鸿沟

在统计学领域，统计学家经常不考虑能得到多少，而是考虑损失了多少，咱们认为损失的负数就是得到。如何衡量损失是一件值得深刻思考的事情。

0x1：例子：只考虑预测结果不考虑不一样预测结果背后的实际损失带来的问题

考虑以下例子：

假设一个气象学家正在预测飓风袭击他所在城市的几率是多少。他估计飓风不袭击该城市的可能性为 99%-100%，且这一律率的可信度为 95%。气象学家很是满意他的计算精度，并建议你们不用疏散。不幸的是飓风袭击了城市，城市被淹没了。

这个典型的例子说明了单纯使用精度去度量结果时存在的缺陷，很显然，这是咱们不少开发朋友包括笔者本身平时在项目中无时不刻遇到的状况。

咱们很幸苦收集了大量打标样本，设计了一套特征工程以及算法模型后，又通过了好久的训练获得了一个预测模型，这时候，预测模型能够获得几千到几万不等的预测结果，其中置信度几率从0.8 ~ 0.99不等，这个时候问题来了，咱们要拿哪些结果进行输出，例如输出给下游业务方或者输出给客户？

是 >= 0.99？那么是否意味着漏报了不少数据？

是 >= 0.9？那么是否意味着产生不少误报？

退一步说， >= 0.99 就必定是100%没误报的吗？咱们设定的这个所谓的 D（>= 0.99）的决策函数，100%不会出问题吗？

形成这个问题的根本缘由是什么呢？

笔者认为根本缘由是：咱们对损失函数的定义还不够完善，如本篇blog的tilte所说，损失函数的本质是”链接统计推断和问题域的桥梁“，因此损失函数就不能单纯只考虑统计推断的几率结果，还要考虑问题域的实际状况，即不一样的决策结果带来的实际现实世界的损失，损失函数必须是承上启下的。

0x2：大体的正确比精确的错误更好

从贝叶斯推断的世界观来看，咱们不能太过于强调精确度的度量，尽管它每每是一个有吸引力和客观的度量，但会忽视了执行这项统计推断的初衷，那就是：推断的结果。

此外，咱们但愿有一种强调“决策收益的重要性”的方法，而不只仅是估计的精确度。

0x3：以寻找最优结果为目标的损失评价函数

在机器学习项目中，咱们更在乎的是获得最佳的预测结果，而不只仅是在全部可能的参数下达到最佳精度。
寻找贝叶斯行动等价于寻找一些参数，这些参数优化的不是参数精度，而是任意某种表现。不过咱们须要先定义表现（损失函数、AUC、ROC、准确率 / 召回率）

Relevant Link:

https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

 

2. 损失函数

在开始讨论这章以前，笔者但愿先抛出一句话：损失函数是链接统计推断和问题所在领域的桥梁。咱们整篇文章会围绕这个观点展开讨论。

损失函数不只仅局限于咱们在机器学习和深度学习书籍里看到的最小二乘损失、交叉熵损失、0-1损失等等。实际上，损失函数是能够根据实际状况，在具体项目中进行自定义的定义，它的核心做用是将现实世界里的问题经过数学方式数字化表达，让损失这个抽象的逻辑概念能够数值化并经过优化算法进行运算和优化。

0x1：损失函数的定义

咱们先来学习统计学和决策理论中的损失函数。损失函数是一个关于真实参数以及对该参数的估计的函数：

损失函数的重要性在于，他们可以衡量咱们的估计的好坏：损失越大，那么损失函数获得的估计就越差。

一个简单而广泛的例子是平方偏差损失函数，这是一种典型的，与偏差的平方成正比的损失函数。在线性回归、无误差统计量计算、机器学习的许多领域都有普遍的应

0x2：常见的损失函数 

1. 0-1损失函数

0-1损失很是直观和简单，非对即错，是一种阶跃函数，在实际项目不多使用到。

2. 平方损失函数 - 非线性损失评价函数

平方损失偏差损失的缺点在于它过于强调大的异常值。这是由于随着估计值的偏离，损失是平方增长的，而不是线性增长的。

对3个单位偏离的惩罚要小于对5个单位误差的惩罚；可是对3个单位偏离的惩罚却不比对1个单位偏离的惩罚大不少，虽然在一些场景下，这两种状况下偏离的差值是相同的。

均方损失的这种不线性，意味着较大的偏差会致使糟糕的结果。在实际项目中，是否可使用须要数据科学家本身根据经验去权衡。

若是咱们的问题是回归分析的场景，咱们是但愿拟合模型能尽可能地归纳数据的分布，离得越远就意味着回归拟合能力越差，由于惩罚也越大，所以可能均方损失的这种非线性特性反而是合理的。

所以要选择损失函数的时候，须要咱们对问题场景要很是的了解，知道选择不一样的损失函数对最终得到的结果意味着什么。

3. 对数损失函数（交叉熵损失函数）

交叉熵损失函数基于信息论中的熵值概念，将损失当作是当前系统不肯定性的度量，预测值和目标值差距越大，意味着不肯定性越大，熵值也就越大，损失也就越大。

4. 绝对值损失函数 - 强调偏离度的线性损失评价函数

一种更稳健的损失函数是偏差的绝对值函数，公式以下：

能够看到，绝对值损失函数不考虑预测值和目标值是”左偏“仍是”右偏“，它只是线性地衡量”预测值和目标值的距离“。

绝对值损失函数是机器学习和稳定统计中常用的损失函数。

笔者插入：能够很容易地看到，绝对损失函数没有对偏离的正负进行区分。but！这个重要吗？截止目前为止，这个问题笔者很难做出回答，我本身也在摸索中，从笔者目前为止为数很少的项目经验来看，两种状况都会遇到，你须要本身灵活判断。

假如你须要的预测的一台机器的恶意网络发包行为的可疑程度，咱们对网络established行为进行建模，进行回归分析，回归模型拟合的值实际上是一个阈值，低于或者高于这个回归值的行为模式都被认为是有必定可疑程度的。可是很显然，从安全业务场景里去定义这个问题时，低于这个阈值被认为是可疑程度相对较低的，高于这个阈值的被认为是可疑程度相对较高的，这个时候，绝对值损失函数就不能100%客观地表征我对这个业务场景的理解了，我可能须要将绝对值损失函数进行改造，改成分段函数，并在负分段和正分段分别赋予不一样的惩罚因子。

咱们在下个小节会讨论的分段损失函数，能够缓解这个问题。

0x3：面向具体业务场景的自定义损失函数

从笔者的理解上来看，目前数据科学家和数据工程师们选择损失函数主要基于以下两方面：

1. 数学计算上的便利性，因为计算机在数学上的便利性，咱们能够自由地设计本身的损失函数；
2. 应用上的稳健性，即损失函数能够客观的对损失进行度量。损失函数确实是客观的，它们都是估计值和真值之差的函数，无论该偏差是正仍是负，和估计的收益对象无关；
3. 损失函数是”可优化的“，这涉及到凸优化理论，理论上说，只有凸函数能够经过计算进行优化，可是值得注意的是，如今对非凸函数的优化也有大量的研究，有兴趣的朋友能够google相关资料；

咱们文章开头说过，损失函数的定义须要承上启下，既要可以对几率问题进行数值化，也要考虑实际问题域的状况。这就要求咱们对常规损失函数进行定制。

再次考虑以上飓风的例子，统计学家也能够这样预测：

飓风袭击的几率在 0%到 1%之间，可是若是发生飓风的损失是巨大的。99%的几率没有飓风，没有飓风带来的收益有限。基于这种损失评价体系，统计学家最后给出的建议可能大相径庭。

很显然，这种改动意味着咱们把关注的重心从更精确的几率估计转到几率估计带来的现实结果上来，即损失函数要更加关注决策结果和业务收益。

1. 分段区间非对称平方偏差损失函数

设想一种场景，预测值低于目标值形成的损失，小于预测值高于目标值。咱们能够用一种非对称的均方偏差损失函数来描述这种关系：

上式代表，这种损失函数对预测值小于目标值的状况，惩罚更严重。在这种损失函数的驱动下，模型会尽可能避免让预测值小于目标值。

这种损失函数形式更适合于估计下个月的网络流量，由于为了不服务器资源分配不足状况，在对服务器资源分配进行估计时要多分配一些。

从这个例子能够看到，损失函数并不必定都是简单的说去衡量模型的预测值和目标值之间的数值差别，而是更加贴近咱们的业务目标。损失函数衡量的是咱们的模型预测结果对咱们的业务目标的贡献程度。

2. 强调参数估计结果更接近 0 或者 1 的损失函数

该损失函数更倾向于估计靠近0或者1位置的参数 𝜃。或者理解为当当目标值为靠近0或者1时，对预测值的一点点偏离都会产生比较大的惩罚。

由于若是真值𝜃靠近 0 或者 1 时，损失函数将变的很是大，除非估计值𝜃̂也 接近 0 或者 1。

基于这种损失函数进行参数优化，模型获得的参数是倾向于靠近 0 或者 1 的。

这种损失函数更适合政治评论家使用，由于他们更倾向于给出 “是/否”的答案（即目标值为0或者1）。

3. 倾向于偏差小的估计的损失函数

该损失函数的取值范围为 0 到 1，这种损失函数代表其使用者不太关心偏差大的估计结果，即预测值和目标值偏离程度并不会显著致使损失增长，这和 0-1 损失函数比较相似，可是在真值附件的惩罚程度又没有 0-1 损失函数那么强烈。0-1损失的惩罚太阶跃了。

4. 带有对某个结果倾向性的损失函数

天气预报人员使用的损失函数经常是带有明显的倾向性的。气象预报员有很强的动力尽可能准确地预报下雨的几率，也一样有动力去错误地暗示可能有雨。为何会出现这种倾向性？

人们更喜欢有所准备，即便可能不会下雨，也不喜欢下雨形成的措手不及。出于这个缘由，预报员倾向于人为地增长降雨几率和夸大这一估计，由于这能带来更大“收益”。

这个和咱们在开头讨论的洪水预测有相似的意思，正是中国那句老话：宁肯信其有，不可信其无。

 

3. 贝叶斯推断理论体系下的损失函数

0x1：样本集数据驱动下的经验损失函数

须要明白的是，咱们以前在讨论损失函数时，都假设咱们知道了真值（目标值），可是实际状况下，损失函数理论下的参数的真实值基本上是不可能已知的，不然，若是咱们已经知道真实值，也就没有必要再去估计了。讨论和真值之间的损失只是为咱们提供一个理论的上界，做为一个benchmark参考系。

所以无论是在传统机器学习/深度学习中，仍是贝叶斯推断中，损失函数都不是针对真是参数值的，而是针对训练样本集的，即经验损失函数。

在这一点上，贝叶斯统计推断和机器学习是同样的，本质上都是训练样本数据驱动的。

0x2：贝叶斯后验分布估计与贝叶斯点估计 - 从几率分布和最大后验估计的角度看后验分布结果

在贝叶斯推断中，把未知参数当作是与先验和后验分布的随机变量。贝叶斯推断不会给出一个具体的估计值，而是该出一个估计值的后验分布。

就后验分布而言，从后验分布获得的样本均可能是真实的参数值（真实的参数值可能就是分布中的某个点）。

而贝叶斯点估计是什么呢？咱们能够理解为在贝叶斯后验估计的几率分布空间中取一个点做为输出结果。这么讲有点抽象，举一个例子说明。

一维高斯函数：

 

上图中的曲线能够理解为对均值和方差【u，σ】的贝叶斯后验估计，它是一个后验几率分布。

而咱们根据最大后验似然估计的原则，选择【u，σ】做为这个分布的最终输出，获得的【u，σ】就是一个贝叶斯估计点。

0x3：贝叶斯估计的指望损失 - 贝叶斯估计损失的理论值

当咱们知道了未知参数的后验分布，就能够基于样本集，计算某种参数估计相关的损失函数了。

可是问题是后验分布的损失函数仍是一个分布，咱们没法作出决策，咱们须要一个肯定的结果或者实值，例如在不少项目中，当咱们经过ML获得一个预测几率值的时候，这时咱们面临一种不肯定性时，咱们每每会经过决策函数（例如一个经验阈值）来将不肯定性提炼成一个单独的肯定性决策动做（0/1）。

咱们须要提炼咱们的后验分布为单一的值（或向量（在多变量的状况下））。若是后验分布的值取得合理，能够避免频率学派没法提供不肯定性的缺陷，同时这种结果的信息更丰富。

回到贝叶斯估计获得的后验估计话题，比起贝叶斯后验估计，咱们更关心的是在某种估计下的指望损失，即一个肯定的实值。

这么作的缘由有不少，首先指望损失有利于解释贝叶斯点估计。现实世界中的系统和机器没法把先验分布做为输入参数。实际上最终关心的是估计结果，后验分布只是计算估计结果的中间步骤， 若是直接给出后验分布是没有太大做用的。

与之相似，咱们须要使后验分布变得陡峭，理想状况是变成一个点。若是后验分布的值取得合理，能够避免频率学派没法提供不肯定性的缺陷，这种结果的信息更丰富。

从贝叶斯后验中选择点， 就是贝叶斯点估计，也即求指望损失。

假设 𝑃(𝜃|𝑋) 是在观测数据 X 的条件下 𝜃 的后验几率，𝜃 的指望损失  为

从后验分布获得 𝑁 个样本 𝜃𝑖 , 𝑖 = 1, ⋯ , 𝑁， 损失函数为𝐿，能够经过大数定理，利用 𝜃̂ 的估计值近似计算出指望损失

笔者插入：这里要注意，指望损失计算的对象贝叶斯估计的理论值，在实际工程中，咱们是没有上帝视角，咱们只有训练样本，所以基于训练样本获得的贝叶斯后验分布以及贝叶斯估计结果，自己就包含了必定的偏置（bias），咱们将基于样本获得的贝叶斯估计结果成为经验损失。

0x2：贝叶斯点估计相比于频率派估计的优点

注意，与 MAP 相比，计算损失函数的指望值利用了后验分布中的更多信息，由于在 MAP 中仅仅是寻找后验分布的最大值，而忽略了后验分布的形状。

MAP 中忽略掉的信息可能使你暴露在长尾部分的风险当中，像飓风这种可能性很小可是存在的风险却很巨大，MAP 将致使估计结果无视参数的未知性。

相似地，频率学派的目标只是最小化化偏差，不考虑与偏差结果相关的损失。频率派方法几乎能够保证永远不会绝对地准确（咱们都是预测模型几乎不多能获得100%的预测结果）。

贝叶斯点估计方法是经过提早计划解决这个问题：若是你的估计将是错误的，那还不如以正确的姿态犯错 --- 模糊地正确性赛过精确的错误。

 

4. 贝叶斯统计实例：优化“价格竞猜”游戏的展品出价

这个小节咱们使用贝叶斯推断中的核心思想：使用基于面向最终结果收益的损失函数，进行策略的优化。

0x1：业务问题场景

咱们的目标是在一个价格竞猜的游戏中取得最大收益，“价格竞猜”的游戏规则以下：

1. 比赛双方争夺竞猜站台上的商品的价格。
2. 每位参赛者看到的商品都不同，价格都是独一无二的，因此并不存在价格各项干扰的状况。
3. 观看后，每位参赛者被要求给出对于本身那套奖品的投标价格。
4. 若是投标价格超过实际价格，投标者被取消获奖资格。
5. 若是投标价格低于真正的价格，且差距在250$之内，投标者得到两套奖品。

游戏的难度在于如何平衡价格的不肯定性，保持你的出价足够低，以便不超过实际价格，同时又不能过低，须要尽可能接近真实价格。

0x2：对问题进行数学建模

很明显，这个问题是一个回归预测问题，解决这个问题的方法有不少种，咱们可使用机器学习回归分析，也能够构建一个深度学习网络进行有监督学习预测，同时，咱们也可使用本章介绍的主题：贝叶斯几率，将咱们的先验知识、训练样本、后验估计结果所有几率化，在几率分布的框架内对问题进行建模分析。

1. 肯定随机变量

无论是机器学习/深度学习/几率统计方法，在建模的时候都暗含了一个要求，即确认随机变量。

因此接下来的问题是，在此次编程建模中，随机变量是什么？

1. 基于历史价格对对待预测价格的先验估计：是输入x；
2. 后验估计分布：是预测值（y）；
3. 对其的信心分布参数：是函数模型的权重参数；
4. 在训练中观察到两个奖项更新后的真实价格：是Y目标值；

2. 肯定先验分布

先验分布的选择，是一个很是复杂的话题，在PAC可学习理论中，先验假设是对假设搜索空间的一种约束，先验越强，这种约束就越强。先验分布的好处是使搜索结果更快且更靠近先验约束，缺点就是若是先验引入的很差，可能会致使最后的搜索结果欠拟合。

这里咱们假设咱们记录了以前的“价格竞猜”比赛，咱们将历史的竞猜结果做为真实价格的先验分布。

咱们假设它遵循一个正态分布：真实价格 ~ Normal（𝜇p, 𝜎p）

经过以前的记录，获得的先验假设𝜇p = 35000，𝜎p = 7500。

笔者思考_1：注意，这里咱们的先验不是一个肯定的实值，而是一个几率分布，在贝叶斯几率编程体系中，读者朋友必定要深入去理解这种”一切皆有可能“的思想，不管是先验知识，仍是后验估计，全部的东西都是一个几率分布，世界不是肯定的，而是几率的。

笔者思考_2：这个问题是一个凭空预测价格的场景，并非提供商品的某些特征而后基于这些特征来预测价格，因此这里其实暗含了一个前提，全部待预测价格的商品，价格都是接近的，收敛在某个价格区间中，因此不管是先验分布仍是信念分布，训练过程本质上都是在尽可能地拟合（靠近）那个价格区间，若是没有这个前提，这个问题是无解的，或者至少不能用简单的一个信念函数来解，请读者朋友仔细体察这点。

3. 决策函数的选择 - 函数模型F() - 信念函数

咱们把先验分布理解为机器学习中的输入x，决策函数至关于F(x)的形式，决策函数负责将先验分布转化为后验几率估计。同时它也是决定咱们应该怎么玩游戏的策略判断依据。

咱们如今有一个关于价格的大概想法（即先验分布估计），但这种猜想可能与真实价格显著不一样（在舞台上压力会成倍增长，你能够看到为何有些出价这么疯狂）。因此咱们须要定义一个信念函数，这个信念函数决定咱们有多大可能性会相信本身的先验判断。

咱们假设咱们关于奖品价格的信念也是正态分布：

Prize i~ Normal（𝜇𝑖, 𝜎𝑖），i = 1，2。这里 i = 1，2 表示每一个参赛者本身的那套奖品只有两个奖品（但这能够扩展到任何数量）。

该套奖品的价格，能够由 Prize1 + Prize2 + w 决定，其中 w 是一个偏差项。

训练的结果是获得一个最匹配训练样本（历史价格分布）的模型参数，即信念函数，这个信念函数是一个连续正态函数，根据不一样的价格给出不一样的信念，基本上来讲，越偏离信念均值中心的价格，信念就越低。

4. 肯定后验估计函数 - 对真实价格估计的模型

首先明确一点，贝叶斯估计是面向几率分布的一种算法，算法获得的后验估计不是一个肯定的实值，而是一个分布。

咱们选择高斯分布做为后验估计的模型，理由是高斯分布很是适合进行决策。高斯分布的均值中心自然就包含了最佳后验几率估计的特性。

同时由于奖品有2个，全部最终的结果是两个分布累加的结果，即后验估计模型函数F(x)是一个复合函数，F(x) = F(x1) + F(x2)。 

0x3：基于PyMC预测奖品价格实例

1. 先验价格分布

让咱们取一些具体的值，假设套装有两个奖品：一趟奇妙的加拿大多伦多之旅，一个可爱的新吹雪机。

咱们对这些奖品的真实价格有一些猜想（先验假设），可是咱们也很是不肯定。按照上面所讨论的，咱们能够用正态分布来表达这种不肯定性。

即吹雪机 ~ Normal（3000，500），多伦多之旅 ~ Normal（12000，3000）

这个先验假设包含了一个事实，咱们相信前往多伦多旅行的真实价格为 12000 美圆，但有68.2% 的几率价格会降低1个标准差，也就是说，有 68.2% 的几率行程价格在 [9000,15000]区间中。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt


norm_pdf = stats.norm.pdf

plt.subplot(311)
x = np.linspace(0, 60000, 200)
# 历史价格先验分布
sp1 = plt.fill_between(x, 0, norm_pdf(x, 35000, 7500),
                       color="#348ABD", lw=3, alpha=0.6,
                       label="historical total prices")
p1 = plt.Rectangle((0, 0), 1, 1, fc=sp1.get_facecolor()[0])
plt.legend([p1], [sp1.get_label()])

plt.subplot(312)
x = np.linspace(0, 10000, 200)
# 吹雪机价格先验分布
sp2 = plt.fill_between(x, 0, norm_pdf(x, 3000, 500),
                       color="#A60628", lw=3, alpha=0.6,
                       label="snowblower price guess")

p2 = plt.Rectangle((0, 0), 1, 1, fc=sp2.get_facecolor()[0])
plt.legend([p2], [sp2.get_label()])

plt.subplot(313)
x = np.linspace(0, 25000, 200)
# 多伦多之旅价格先验分布
sp3 = plt.fill_between(x, 0, norm_pdf(x, 12000, 3000),
                       color="#7A68A6", lw=3, alpha=0.6,
                       label="Trip price guess")
plt.autoscale(tight=True)
p3 = plt.Rectangle((0, 0), 1, 1, fc=sp3.get_facecolor()[0])
plt.legend([p3], [sp3.get_label()])

plt.show()

2. PYMC训练及预测后验几率估计

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm

norm_pdf = stats.norm.pdf

plt.subplot(311)
x = np.linspace(0, 60000, 200)
# 历史价格先验分布
sp1 = plt.fill_between(x, 0, norm_pdf(x, 35000, 7500),
                       color="#348ABD", lw=3, alpha=0.6,
                       label="historical total prices")
p1 = plt.Rectangle((0, 0), 1, 1, fc=sp1.get_facecolor()[0])
plt.legend([p1], [sp1.get_label()])

plt.subplot(312)
x = np.linspace(0, 10000, 200)
# 吹雪机价格先验分布
sp2 = plt.fill_between(x, 0, norm_pdf(x, 3000, 500),
                       color="#A60628", lw=3, alpha=0.6,
                       label="snowblower price guess")

p2 = plt.Rectangle((0, 0), 1, 1, fc=sp2.get_facecolor()[0])
plt.legend([p2], [sp2.get_label()])

plt.subplot(313)
x = np.linspace(0, 25000, 200)
# 多伦多之旅价格先验分布
sp3 = plt.fill_between(x, 0, norm_pdf(x, 12000, 3000),
                       color="#7A68A6", lw=3, alpha=0.6,
                       label="Trip price guess")
plt.autoscale(tight=True)
p3 = plt.Rectangle((0, 0), 1, 1, fc=sp3.get_facecolor()[0])
plt.legend([p3], [sp3.get_label()])

# plt.show()

# 吹雪机和多伦多之旅的价格猜想先验分布
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 对价格的信念分布就是历史价格的先验分布
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)


mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


# predict
x = np.linspace(5000, 40000)
plt.plot(x, stats.norm.pdf(x, 35000, 7500), c="k", lw=2,
         label="prior dist. of suite price")

_hist = plt.hist(price_trace, bins=35, normed=True, histtype="stepfilled")
plt.title("Posterior of the true price estimate")
plt.vlines(mu_prior, 0, 1.1 * np.max(_hist[0]), label="prior's mean",
           linestyles="--")
plt.vlines(price_trace.mean(), 0, 1.1 * np.max(_hist[0]),
           label="posterior's mean", linestyles="-.")
plt.legend(loc="upper left")

plt.show()

从上图能够看到，基于吹雪机和多伦多旅行的价格先验，整体价格从历史先验的35000降低到了20000（图中蓝色色块）。

3. 结合实际业务场景的损失函数

从最大后验的角度来看，咱们此时已经解决问题了，直接给出最终出价便可。可是在贝叶斯统计领域，咱们有更多的信息，咱们应该将损失归入咱们的最终出价中。
咱们须要使用结合业务场景的损失函数来找到最佳出价，即对于咱们的损失函数来讲是最优解。

def showcase_loss(guess, true_price, risk=80000):
    if true_price < guess:
        return risk
    elif abs(true_price - guess) <= 250:
        return -2*np.abs(true_price)
    else:
        return np.abs(true_price - guess - 250)

其中 riks 是一个参数，代表若是你的猜想高于真正的价格的损失程度。
这里选择了一个值: 80000，风险较低意味着你更能忍受出价高于真实价格的状况。
若是咱们出价低于真实价格，而且差别小于250，咱们将得到两套奖品；不然，当咱们出价比真实价格低，咱们要尽量接近，所以其损失是猜想和真实价格之间距离的递增函数。

接下来的问题是，这个 risk 参数如何选择呢？

对于每一个可能的出价，咱们计算与该出价相关联的指望损失，经过改变 risk 参数来看它如何影响咱们的最终损失。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm

norm_pdf = stats.norm.pdf

# plt.show()

# 吹雪机和多伦多之旅的价格猜想先验分布
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 对价格的信念分布就是历史价格的先验分布
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)

mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


def showdown_loss(guess, true_price, risk=80000):
        loss = np.zeros_like(true_price)
        ix = true_price < guess
        loss[~ix] = np.abs(guess - true_price[~ix])
        close_mask = [abs(true_price - guess) <= 250]
        loss[close_mask] = -2 * true_price[close_mask]
        loss[ix] = risk
        return loss


guesses = np.linspace(5000, 50000, 70)
risks = np.linspace(30000, 150000, 6)
expected_loss = lambda guess, risk: \
    showdown_loss(guess, price_trace, risk).mean()

for _p in risks:
    results = [expected_loss(_g, _p) for _g in guesses]
    plt.plot(guesses, results, label="%d" % _p)

plt.title("Expected loss of different guesses, \nvarious risk-levels of \
overestimating")
plt.legend(loc="upper left", title="Risk parameter")
plt.xlabel("price bid")
plt.ylabel("expected loss")
plt.xlim(5000, 30000)

plt.show()

 

最大限度地减小咱们的损失是咱们的目标，这对应于上面图中的每条曲线的最小值点。更正式地说，咱们但愿经过寻找如下公式的解来减小咱们的指望损失

指望损失的最小值被称为贝叶斯行动。咱们能够基于SciPy的fmin api进行智能搜索，寻找任一单变量或多变量函数的极值（不必定是全局极值）

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop

norm_pdf = stats.norm.pdf

# plt.show()

# 吹雪机和多伦多之旅的价格猜想先验分布
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 对价格的信念分布就是历史价格的先验分布
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)

mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


ax = plt.subplot(111)


def showdown_loss(guess, true_price, risk=80000):
        loss = np.zeros_like(true_price)
        ix = true_price < guess
        loss[~ix] = np.abs(guess - true_price[~ix])
        close_mask = [abs(true_price - guess) <= 250]
        loss[close_mask] = -2 * true_price[close_mask]
        loss[ix] = risk
        return loss

guesses = np.linspace(5000, 50000, 70)
risks = np.linspace(30000, 150000, 6)
expected_loss = lambda guess, risk: \
    showdown_loss(guess, price_trace, risk).mean()


for _p in risks:
    _color = next(ax._get_lines.prop_cycler)
    _min_results = sop.fmin(expected_loss, 15000, args=(_p,),disp = False)
    _results = [expected_loss(_g, _p) for _g in guesses]
    plt.plot(guesses, _results , color = _color['color'])
    plt.scatter(_min_results, 0, s = 60, \
                color= _color['color'], label = "%d"%_p)
    plt.vlines(_min_results, 0, 120000, color = _color['color'], linestyles="--")
    print("minimum at risk %d: %.2f" % (_p, _min_results))

plt.title("Expected loss & Bayes actions of different guesses, \n \
various risk-levels of overestimating")
plt.legend(loc="upper left", scatterpoints=1, title="Bayes action at risk:")
plt.xlabel("price guess")
plt.ylabel("expected loss")
plt.xlim(7000, 30000)
plt.ylim(-1000, 80000)

plt.show()

当咱们下降风险阈值（不那么担心出价太高），咱们能够提升咱们的出价，想要更接近真实价格。

这个例子中，由于维度较低，咱们能够肉眼找到极值点，可是在高维函数中，肉眼是没法找到的极值的。
这里列举几个损失函数的贝叶斯行动能够用显式公式表达的例子：

1. 使用均方损失: 贝叶斯行动是后验分布的均值；
2. 当后验分布的中位数将绝对指望损失函数最小化时，用样本的中位数来接近是很是准确的；
3. MAP估计是某个损失函数收缩到0-1损失的解；

Relevant Link:

https://www.zhihu.com/question/21134457
http://www.xuyankun.cn/2017/05/13/bayes/
http://mindhacks.cn/2008/09/21/the-magical-bayesian-method/

 

5. 贝叶斯行动预测的另外一个例子 - 金融预测

0x1：问题场景

咱们须要设计一个模型，对某只股票将来的价格进行预测，咱们的利润和损失将直接依赖于基于预测价格的行动。
很显然，若是预测将来价格是上行的(价格上涨)，而实际也上涨，那不管是精确预测仍是预测结果高于实际价格，至少是盈利的；可是若是股票价格实际下行了，预测倒是上行，那么不管预测结果是靠近仍是偏离，都是亏损的，区别只是亏多亏少而已，固然，亏少确定比亏多要好。
同时还有一个约束条件，做为金融机构来讲，稳健的收益回报比大比例的盈利和亏损更被看重，因此即便是预测上行对了，也尽可能会避免过大的偏离实际值。

0x2：模型设计

1. 肯定损失评价函数

如何衡量模型的预测结果和将来实际价格之间差距的损失？假如咱们使用平方偏差损失函数，那么对正负号是不加以区分的，预测值 -0.01 和 0.03 和实际价格 0.02 之间的惩罚相同: 

(0.01 - (-0.01))2 = (0.01 - 0.03)2 = 0.0004

若是咱们基于采用了这种损失策略的模型的预测下了堵住，那么 0.03 的预测值会使你赚钱，而 -0.01 的预测值会使你赔钱，但损失函数没法提示这一点。咱们须要一个更好的损失函数，即考虑了预测价格的正负号和真正的价值的损失函数。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop


def stock_loss(true_return, yhat, alpha=100.):
    if true_return * yhat < 0:
        # opposite signs, not good
        return alpha * yhat ** 2 - np.sign(true_return) * yhat \
            + abs(true_return)
    else:
        return abs(true_return - yhat)


true_value = .05
pred = np.linspace(-.04, .12, 75)

plt.plot(pred, [stock_loss(true_value, _p) for _p in pred],
         label="Loss associated with\n prediction if true value = 0.05", lw=3)
plt.vlines(0, 0, .25, linestyles="--")

plt.xlabel("prediction")
plt.ylabel("loss")
plt.xlim(-0.04, .12)
plt.ylim(0, 0.25)

true_value = -.02
plt.plot(pred, [stock_loss(true_value, _p) for _p in pred], alpha=0.6,
         label="Loss associated with\n prediction if true value = -0.02", lw=3)
plt.legend()
plt.title("Stock returns loss if true value = 0.05, -0.02")

plt.show()

上图中咱们能够看到几个关键的特征：

1. 损失函数反应了用户并不想猜想符号，猜错符号的损失远大于猜对符号；
2. 无论是否猜对符号，用户都不想大幅度猜错，即大幅度偏离真实值。这也体现了金融机构应对下行风险(预测方向是错的，量级很大)和上线风险(预测方向是正确的，量级很大)的态度是类似的。二者都被视为危险的行为而不被孤立。所以，当咱们进一步远离真实价格，咱们的损失会增长，但在正确方向上的极端损失相对错误方向上的较小。

2. 肯定函数模型F()

咱们将会在被认为可以预测将来回报的交易信号上作一个回归。数据是人工虚构的，由于绝大多数时候金融数据都不是线性的。

下图中，咱们沿着最小方差线画出了数据分布状况：

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop


# Code to create artificial data
N = 100
X = 0.025 * np.random.randn(N)
Y = 0.5 * X + 0.01 * np.random.randn(N)

ls_coef_ = np.cov(X, Y)[0, 1] / np.var(X)
ls_intercept = Y.mean() - ls_coef_ * X.mean()

plt.scatter(X, Y, c="k")
plt.xlabel("trading signal")
plt.ylabel("returns")
plt.title("Empirical returns vs trading signal")
plt.plot(X, ls_coef_ * X + ls_intercept, label="Least-squares line")
plt.xlim(X.min(), X.max())
plt.ylim(Y.min(), Y.max())
plt.legend(loc="upper left")

plt.show()

线性回归的函数形式以下：

对一个特定的交易信号 x，回报的分布有以下形式：

对于给定的损失函数，咱们但愿找到：

下图中咱们对比了在平方损失函数和咱们自定义损失函数条件下，不一样的交易信号的贝叶斯行动。

从上图中能够看出：

1. 当交易信号接近为0时，正负回报均可能出现，咱们最好的行为是预测接近为0，也就是说处于中立；
2. 只有当咱们都很是自信时，咱们才进场下注，当对不肯定性感到不安时，选择不做为；
3. 当信号变得愈来愈极端时，咱们对正负回报的预测愈来愈自信，自定义损失函数的模型将会收敛到最小二乘；

具备以上特色的模型被称为系数预测模型。稀疏模型不是想要想方设法地去拟合数据，稀疏模型试图找打针对咱们定义的股票损失的最好预测，而最小二乘只试图找到相对于平方损失下的数据的最佳拟合。

Relevant Link: