TowardsDataScience 博客中文翻译 2020（九百五十六）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

VAR 和面板数据模型——多元预测技术的发电站

原文：https://towardsdatascience.com/var-and-panel-data-models-the-powerhouse-of-multivariate-forecasting-techniques-22b8d8888141?source=collection_archive---------28-----------------------

艾伦蔡在 Unsplash 上的照片

时间序列

超越单变量时间序列预测

我们来做一个思维实验。

你是一家专门从事销售预测的预测分析公司。您的客户是一家大型食品连锁店，他要求您预测未来 30 天的销售量。

现在，您将寻找哪些信息来进行预测？

当然是过去的销售数据。假设你得到了他们过去两年的销售数据。所以你的预测模型变成了这样:

未来销售额= f(过去销售额)

现在想想，你能有多自信地说，只用这一个历史数据集，你就能做出一个好的、合理的预测？

在现实世界中当然不是这样。在现实世界中，预测历史销售额只是一条信息，但你还会寻找商店位置、法定假日、任何特殊事件、广告支出或促销活动等，事实上，每一条影响销售的信息。

这些信息的每一部分实质上都成为了预测模型中的一个变量:

未来销售额= f(过去的销售额、特殊事件、假期、促销……)

这就是多元预测的力量所在。现实世界是相互联系的，一个过程的结果通常是多个输入和动作的结果。因此，难怪数据科学家会寻找与构建预测模型相关的每一条信息。

在本文中，我将讨论两个强大的多变量时间序列预测模型——向量自回归(VAR)和面板数据模型——并用两种不同的编程语言——Python 和 r——的代码片段演示它们的应用。

所以让我们开始吧。

一.向量自回归模型

a)什么是向量自回归？

在定义向量自回归之前，让我们了解一下什么是自回归。顾名思义，它是一个变量对自身的回归，对它过去的值的回归——就像我刚刚在引言中写的那样。在这种情况下，独立变量是被预测的同一个单变量数据序列的过去值。

向量自回归概括了这一单变量概念，并允许在模型中包含额外的相关变量。在这种情况下，每个变量都使用其自身的过去值(滞后值)以及外部因素的滞后值进行预测。

a)风险值模型的工作原理

在单变量预测技术中，未来值被预测为过去值的函数。所以，如果你想预测一个城市在 2021 年的人口，你需要一个简单的单变量预测方法:

但如果你想在更长的时间范围内做出预测(比如 2050 年的人口)，我们知道还有许多其他因素需要考虑——出生率、死亡率、教育、收入等等。这些因素被称为外源因素或协变量。在短期内，这些外部因素的影响不是很大，所以我们可以忽略它们，但从长期来看，它们的影响是复合的，并产生了可测量的差异。

VAR 就是这样一种在建模过程中融入外生变量的多变量预测模型。在此过程中，因变量使用其自身的过去(滞后)值以及外部因素的滞后值进行预测。

c)用 Python 实现

为了演示它是如何工作的，我使用了一个由 14 个不同变量组成的美国宏观经济数据集。这是从 1959 年到 2009 年的季度数据。然而，出于演示的目的，我只使用了两个变量:国内生产总值(realgdp)和个人消费支出(realcons)。您可以按照下面代码片段中的代码访问数据集。

用于 VAR 实施的完整数据集的几行

我将使用一个 Python 环境来利用其非常丰富的statsmodels库，这是一个非常灵活的 VAR 实现选项的来源。演示的数据集也来自这个 API。

所以这个练习的目标是对未来的实际 gdp 和实际成本进行 5 个时间步的预测。

二。面板数据模型

a)什么是面板数据？

如果你之前没有听说过面板数据你并不孤单。但是我确信你以前见过这种数据，你只是不知道它有一个名字。以下是面板数据的正式定义:

一个 面板数据 是一个观察值的多维数据，随着时间的推移被反复测量。

简单地说，一个面板数据是对同一个对象或个体的重复观察。看看下面这些虚构的数据，它们可能比阅读定义更直观:

图:假想面板数据表

b)面板数据模型如何工作？

从概念的角度来看，面板数据很像多元 OLS 回归，但比正常的 OLS 更有效。使用上述数据集，假设您对以下问题感兴趣:

城市支出增加了入学人数吗？

通过运行 OLS，我们可以快速检查城市支出和学校招生之间是否存在关联，并构建一个模型:

支出= f(学校招生)

然而，数据集不是普通的数据框架；这是一个面板数据，有两个额外的列——城市名称和年份——我们也可以在模型中使用。下面是三种不同的方法来构建模型，每种方法都增加了额外的可变效果。

汇集 OLS 模式

固定特效模式

随机效果模型

c)在R中的实施

我喜欢 Python，但这是我在 r 上看到的为数不多的情况之一。它有一个强大的库叫做plm,你可以用它只用几行代码构建面板数据模型。而且非常直观。

下面我正在创建一个关于随机效果模型的演示，但是你可以通过将模型类型(第 16 行)从“随机”改为“池化”或“内化”以及一些其他的建模变化来玩代码。

图:模型输出汇总

最后的想法

本文的目的是展示两种强大的多元预测技术——向量自回归和面板数据模型——在现实应用中的价值。为什么要多元预测？因为单变量预测有其自身的局限性。认为未来完全依赖于过去是一个过于简单的假设。

建立模型的目的是为了重现和模仿现实世界。现实世界是复杂且相互关联的，因此这需要在我们的模型中得到反映。

使用蒙特卡罗模拟的投资组合风险管理

原文：https://towardsdatascience.com/var-calculation-using-monte-carlo-simulations-40b2bb417a67?source=collection_archive---------7-----------------------

如何使用 Python 和蒙特卡罗模拟自动计算风险价值(VaR)来管理投资组合或股权和股票的金融风险

在 Unsplash 上由 Austin Distel 拍摄的照片

金融和投资组合风险管理中的 VaR？

VaR 是“风险价值”的首字母缩写，是许多公司和银行用来确定其公司内部金融风险水平的工具。风险值是为一家公司的投资计算的，或者可能是为了检查由银行财富管理部门或精品公司管理的投资组合的风险水平。

该计算可以被认为是孤立的统计测量。也可以简化为下面的语句示例-

VaR 是在一定概率水平(置信区间)下发生的最小损失或在一定概率水平下实现的最大损失。

图片来源——SP 咨询 LLP 公司

上图显示了一家公司在α % 信心下可能面临的最大损失。在个人层面上，VaR 可以帮助你预测或分析你的投资组合可能面临的最大损失——这是我们将很快分析的内容。

蒙特卡洛模拟

蒙特卡洛模型是斯坦尼斯瓦夫·乌拉姆和约翰·诺伊曼的发明，他们在第二次世界大战后开发了这个模型。该模型是以摩纳哥的一个赌博城市命名的，这是因为赌博中面临的机会和随机遭遇。

蒙特卡洛模拟是一种概率模型，它生成与经济因素(预期回报、波动性——在基金投资组合的情况下)一起使用的随机变量，以预测大范围内的结果。虽然不是最准确的，但该模型通常用于计算风险和不确定性。

我们现在将使用蒙特卡洛模拟来为我们的资产组合生成一组预测回报，这将帮助我们找出我们投资的风险值。

用 Python 计算 VaR

我们将首先通过导入所需的库和函数来设置笔记本

#Importing all required libraries
#Created by Sanket Karve
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pandas_datareader as web
from matplotlib.ticker import FuncFormatter
!pip install PyPortfolioOpt
#Installing the Portfolio Optimzation Libraryfrom pypfopt.efficient_frontier import EfficientFrontier
from pypfopt import risk_models
from pypfopt import expected_returns
from matplotlib.ticker import FuncFormatter

为了我们项目的目的，我考虑了过去两年的‘FAANG’股票。

tickers = ['GOOGL','FB','AAPL','NFLX','AMZN']
thelen = len(tickers)
price_data = []
for ticker in range(thelen): prices = web.DataReader(tickers[ticker], start='2018-06-20', end = '2020-06-20', data_source='yahoo') price_data.append(prices[['Adj Close']])df_stocks = pd.concat(price_data, axis=1)
df_stocks.columns=tickers
df_stocks.tail()

下一步，我们将计算每项资产的投资组合权重。我已经通过使用为达到最大夏普比率而计算的资产权重做到了这一点。我已经发布了下面计算的代码片段。

#Annualized Return
mu = expected_returns.mean_historical_return(df_stocks)
#Sample Variance of Portfolio
Sigma = risk_models.sample_cov(df_stocks)#Max Sharpe Ratio - Tangent to the EF
from pypfopt import objective_functions, base_optimizer
ef = EfficientFrontier(mu, Sigma, weight_bounds=(0,1)) #weight bounds in negative allows shorting of stocks
sharpe_pfolio=ef.max_sharpe() #May use add objective to ensure minimum zero weighting to individual stocks
sharpe_pwt=ef.clean_weights()
print(sharpe_pwt)

最大夏普比率的资产权重

资产权重将用于计算预期投资组合回报。

#VaR Calculation
ticker_rx2 = []#Convert Dictionary to list of asset weights from Max Sharpe Ratio Portfoliosh_wt = list(sharpe_pwt.values())
sh_wt=np.array(sh_wt)

现在，我们将把投资组合的股票价格转换成累计回报，这也可以被认为是该项目的持有期回报(HPR)。

for a in range(thelen): ticker_rx = df_stocks[[tickers[a]]].pct_change()
  ticker_rx = (ticker_rx+1).cumprod()
  ticker_rx2.append(ticker_rx[[tickers[a]]])ticker_final = pd.concat(ticker_rx2,axis=1)
ticker_final

#Plot graph of Cumulative/HPR of all stocksfor i, col in enumerate(ticker_final.columns):
  ticker_final[col].plot()plt.title('Cumulative Returns')
plt.xticks(rotation=80)
plt.legend(ticker_final.columns)#Saving the graph into a JPG file
plt.savefig('CR.png', bbox_inches='tight')

现在，我们将挑选出每项资产的最新 HPR，并用*将收益乘以计算出的资产权重。点()*功能。

#Taking Latest Values of Returnpret = []
pre1 = []
price =[]for x in range(thelen):
  pret.append(ticker_final.iloc[[-1],[x]])
  price.append((df_stocks.iloc[[-1],[x]]))pre1 = pd.concat(pret,axis=1)
pre1 = np.array(pre1)
price = pd.concat(price,axis=1)
varsigma = pre1.std()
ex_rtn=pre1.dot(sh_wt)print('The weighted expected portfolio return for selected time period is'+ str(ex_rtn))#ex_rtn = (ex_rtn)**0.5-(1) #Annualizing the cumulative return (will not affect outcome)price=price.dot(sh_wt) #Calculating weighted value
print(ex_rtn, varsigma,price)

计算出预期投资组合回报和波动性(预期回报的标准差)后，我们将建立并运行蒙特卡洛模拟。我用了 1440(一天中的分钟数)的时间进行了 10，000 次模拟运行。时间步长可以根据要求而改变。我使用了 95%的置信区间。

from scipy.stats import normimport mathTime=1440 #No of days(steps or trading days in this case)lt_price=[]final_res=[]for i in range(10000): #10000 runs of simulation daily_return=                     (np.random.normal(ex_rtn/Time,varsigma/math.sqrt(Time),Time))
  plt.plot(daily_returns)plt.axhline(np.percentile(daily_returns,5), color='r', linestyle='dashed', linewidth=1)plt.axhline(np.percentile(daily_returns,95), color='g', linestyle='dashed', linewidth=1)plt.axhline(np.mean(daily_returns), color='b', linestyle='solid', linewidth=1)plt.show()

一天 1440 分钟的回报范围|红色虚线-最小损失|绿色虚线-最小收益

可视化回报的分布图为我们呈现了下面的图表

plt.hist(daily_returns,bins=15)plt.axvline(np.percentile(daily_returns,5), color='r', linestyle='dashed', linewidth=2)plt.axvline(np.percentile(daily_returns,95), color='r', linestyle='dashed', linewidth=2)plt.show()

打印上限和下限的精确值，并假设我们的投资组合价值为 1000 美元，我们将计算出应该保留以弥补最小损失的资金量的估计值。

print(np.percentile(daily_returns,5),np.percentile(daily_returns,95)) #VaR - Minimum loss of 5.7% at a 5% probability, also a gain can be higher than 15% with a 5 % probabilitypvalue = 1000 #portfolio valueprint('$Amount required to cover minimum losses for one day is ' + str(pvalue* - np.percentile(daily_returns,5)))

每天的最小损失为 1.14%，概率为 5%

得出的金额将意味着弥补你每天损失所需的金额。结果也可以解释为你的投资组合以 5%的概率面临的最小损失。

结论

上述方法显示了我们如何计算投资组合的风险价值(VaR)。复习使用现代投资组合理论(MPT) 计算一定投资金额的投资组合，将有助于巩固您对投资组合分析和优化的理解。最后，VaR 与蒙特卡洛模拟模型相结合，也可用于通过股价预测损益。这可以通过将生成的每日回报值乘以相应股票的最终价格来实现。

所有表达的观点都是我自己的。这篇文章不被视为专家投资建议。

随时联系我LinkedIn|Twitter给我留言

可变大小的视频小批量

原文：https://towardsdatascience.com/variable-sized-video-mini-batching-c4b1a47c043b?source=collection_archive---------38-----------------------

一种处理帧数不等的视频的批处理算法

亨特·哈里特在 Unsplash 上的照片

T 训练和测试高效机器学习模型最重要的一步是收集大量数据并使用这些数据有效训练它的能力。小批量有助于解决这个问题，因为它能够在每次迭代中使用一小部分数据进行训练。但是，随着大量的机器学习任务在视频数据集上执行，出现了有效批处理不等长视频的固有问题。大多数方法依赖于将视频修剪为相等的长度，以便在迭代期间可以提取相同数量的帧并发送给批处理。但在我们需要每一帧的信息来有效预测一些事情的情况下，尤其是在自动驾驶汽车和动作识别的情况下，这不是特别有帮助。这就是创建小批量视频数据集的动机，该数据集可以处理不同长度的视频。

我使用 LoadStreams 作为由Glenn Jocher(https://github.com/ultralytics/yolov3)创建的 Yolov3 的 Pytorch 模型的基础，在此基础上我对我的 LoadStreamsBatch 类进行建模。

类初始化

在 init 函数中，我接受四个参数。虽然 img_size 与原始版本相同，但其他三个参数定义如下:

• sources: 它接受一个目录路径或一个文本文件作为输入。
• **批量大小:**这是所需的最小批量的大小
• subdir_search: 该选项可以切换，以确保当目录作为 sources 参数传递时，在所有子目录中搜索相关文件

我首先检查 sources 参数是一个目录还是一个文本文件。如果是一个目录，我会读取目录中的所有内容(如果 subdir_search 参数为 True ，子目录也包括在内)，否则我会读取文本文件中视频的路径。视频的路径存储在列表中。维护指针 cur_pos 以跟踪列表中的当前位置。

列表以 batch_size 为最大值进行迭代，并进行额外检查以跳过错误视频或不存在的视频。这些被发送到信箱功能，以调整图像的大小，并将所有组件堆叠在一起，这与原始版本没有变化，除非所有视频都有问题/不可用。

函数定期检索帧

由线程调用的 update 函数有一个小的变化，我们额外存储了默认的图像大小，用于所有视频都被拾取进行处理，但由于长度不等，一个在另一个之前完成的情况。当我解释代码的下一部分，也就是 next 函数时，会更清楚。

迭代程序

在获得初始运行的所有帧之后，用于处理的帧可以在每个视频的迭代的不同点结束。如果该帧存在，它像往常一样被传递给信箱函数。在帧为 None 的情况下，这意味着视频已经被完全处理，我们检查列表中的所有视频是否已经被处理。如果有更多视频需要处理，指针 cur_pos 用于获取下一个可用视频的位置。
在列表中不再有视频被拾取，但某些视频仍在处理的情况下，一个空白帧被发送到其他批次组件，即它根据其他批次中的剩余帧动态地调整视频的大小。

结论

由于花费在数据收集和数据预处理上的时间量，我相信这可以通过减少花费在使视频适合模型上的时间而有所帮助，并且我们可以更专注于使模型适合数据。

我在这里附上完整的源代码。希望这有所帮助！

注入差异的思维

原文：https://towardsdatascience.com/variance-infused-thinking-49f3e780890d?source=collection_archive---------52-----------------------

Riho Kroll 在 Unsplash 上拍摄的照片

要成为一名成功的数据科学家或分析从业者，我们需要练习从概率的角度看待世界

历史看起来是确定的，而且经常被这样解释。本世纪初的科技泡沫显然是一个泡沫，正如所有泡沫都会发生的那样，它最终破裂了。20 世纪 80 年代的日本房地产，2008 年的美国房地产，20000 美元的比特币，杂草股等等也是如此。在崩盘后的几年里，历史学家关注的是风险有多明显，最终的内爆有多不可避免。

成功的故事也是如此。迈克尔·乔丹作为篮球运动员的无与伦比的成功基本上是命中注定的，这要归功于他独特的运动能力和竞争力的完美结合。谷歌和亚马逊世界一流的算法和执行能力意味着万亿美元的市值是不可避免的。

我可能有点夸大其词，但是很多人看待历史的一个主要缺陷是他们通过决定论的视角来看待历史。部分原因是历史是线性的——只有一系列事件导致了今天，而这一系列事件就是历史书上所描述的。历史书不会深究如果，比如如果四年前特朗普输给了克林顿会怎样；相反，他们关注(可能过度关注)过去和现在。

决定性思维的陷阱

我经常想知道为什么概率被认为是数据科学家的必修知识。日复一日，我很少使用概率——在我的日常工作中，我更多地使用统计学、矩阵代数和机器学习。我得出的答案是:概率很重要，因为数据科学家需要学习概率思维。概率性思维意味着认识到结果的序列不是确定的——而是每一个结果本身都是因果关系(或相关性)和机会的混合产物。

下表显示了这种差异。假设我们观察到以下事件序列(只观察了一次):A，然后 B，然后 C。观察世界的确定性方法会得出结论，A 导致 B，B 导致 C。在未来，如果我们再次看到 A，我们会假设 C 就在眼前。这种想法有三个陷阱:

1. 我们在不该推断的时候推断出强因果关系(特别是考虑到我们只有一个样本)。
2. 我们没有意识到机遇对结果的巨大影响。
3. 我们只根据结果来判断过程。

**Deterministic:** A ----> B ----> C

另一方面，一种看待事物的概率方法认为，如果 A 出现(也有可能 A 不会出现)，那么它后面可能是 B 或 E(或其他)。无论发生哪种情况，B 或 E，都可能导致一系列后果。因此，即使 A -> B -> C 是我们所观察到的，这并不意味着下一次我们再次看到 A 时，我们应该期待 C。

**Probabilistic:** if A ----> B ----> C
         \       \----> D
          \   
           \--> E ----> F
            \    \----> G
             \    \---> H
              \
               \------> Other

相反，实际发生的事情，C，可能是最不可能的结果。例如，不久前赫兹(一家破产公司)的股票风靡一时。比方说，我买了一些赫兹看涨期权(对股票的杠杆化押注)，几天后卖出，获得了可观的利润。**同样为了完整起见，假设 Hertz 的股票上涨了 100%，我的期权将上涨 300%——因此，如果我下注 1000 美元，Hertz 翻倍，我将净赚 3000 美元。**赚钱后，我的大脑大概会经历以下想法:

Buy HTZ calls ----> Make money ----> I'm a good investor

但这将是一个错误。我甚至可能再次走运，在其他一些蹩脚的破产公司身上下注赚钱；但在某个时候，小鸡会回家休息，我会失去一切。这是因为我没有诚实地思考我的赫兹赌注的可能结果范围，也就是结果锥(我也没有考虑如果我投资其他地方会发生什么)。相反，我过于关注积极的结果。更现实和全面的评估将实现以下目标:

• 赫兹是一家破产的公司。因此，如果提振股票的动物精神突然蒸发，就没有实际的商业价值来保护我的投资不成为零。
• 由于标的股票的波动性，赫兹期权非常昂贵，其每日波动幅度超过 100%(上下波动)。这意味着股票必须大幅上涨，期权才能兑现——这不是一个很可能的结果。

我们可以用概率画出一个更真实更完整的结果树(下图)。既然我们已经描绘出了可能的结果范围，我们可以看到我的赌博是多么愚蠢。我有 5%的机会扭亏为盈，95%的机会亏损(因为你失去了期权到期时你支付的全部溢价)。此外，由于期权非常昂贵(由于波动性)，对我有利的巨大波动的回报只有 3:1——但基于波动的低赔率，我应该要求 20:1 的回报。考虑回报也很重要——即使回报的概率很小，如果回报的大小超过了补偿，那么期望值仍然可以是正的。

所以总而言之，不，我根本没有做好投资，我只是经历了一场幸运的赌博——这绝对是一个运气掩盖了可怕过程的案例，从长远来看，这是不可持续的。

Buy HTZ calls ---> Stock up a lot (Pr. =5%) ----> Options up
          \
           \-----> Stock up a bit (Pr. =15%) ---> Options worthless
            \
             \---> Stock no change (Pr. =20%) --> Options worthless
              \
               \-> Stock down (Pr. =60%) -------> Options worthless

以这种方式陈述时，听起来非常明显——如此明显，你可能会想，谁会愚蠢到这样想呢？但是很多人都这样做。我想到的最明显的例子是某个总统因股市飙升而邀功:

Becomes president ----> Makes economy great ----> Stocks go up

每当他在推特上发布股价创纪录上涨或创历史新高的消息时，就会有成千上万的人陷入确定性思维的陷阱——他们开始相信股市肯定上涨了，因为经济表现良好；多亏了总统，经济发展得很好。但他们没有认识到以下事件的替代顺序和同样可信的解释:

• 股票市场通常与经济脱节。即使在经济低迷的时候，股票有时也会因为纯粹的机会而上涨(股票价格经常以违背理性和解释的方式波动)。因此，总统可能只是运气好(尽管他对经济造成了损害，但股票却不太可能上涨)，他正试图为这种运气邀功。这里有一个更现实的结果树(看起来很明显，我知道)，它也包括了可能发生的情况:

Becomes president ----> Hurts economy ----> Stocks go up anyway
                        with Trade War
                        and pandemic 
                        response                                  
                                  \-------> Stocks go sideways
                                   \
                                    \-----> Stocks go down

• 股票可能会上涨，因为自 2008 年崩盘以来，前所未有的大量流动性被注入金融系统，，这与现任总统无关(与央行政策关系更大)。此外，重要的是要认识到，即使美联储的印钞机全速运转，股价上涨也不是唯一的结果。这可能是最有可能的结果，但也有其他可能的和不太有利的结果。

Fed prints money ----> Stocks go up
                \
                 \---> Stocks go sideways
                  \
                   \-> Stocks go downPresident happens to be elected at the same time (unrelated)

正确的思考方式

从分析角度思考未来的正确方式是问:

如果我把这个游戏玩很多次，可能会有什么样的结果？

把未来想象成一个圆锥形的结果

我们应该总是从分布的角度来思考，把未来想象成一个以我们行动的预期结果为中心的潜在结果的圆锥(或任何最合理地接近分布的形状)。**此外，我们还需要考虑如果不采取行动的后果。**这就是像 A/B 测试这样的受控实验背后的原理——一个实验有一个控制组和一个实验组，适当地隔离出要测试的因素(或行为)。

受控实验类似于模拟，允许我们反复观察如果我们采取行动(上图中的绿色圆锥)和如果我们不采取行动(红色圆锥)会发生什么。如果我们实验的样本量足够大，那么我们就可以对可能的结果范围有一个很好的认识，这使我们能够正确地评估所考虑的行动的预期效果(它是否使我们朝着更理想的平均结果前进？)和围绕行动效果的不确定性()它是否扩大了圆锥并引入了过度的风险？)。

不幸的是，在现实世界中，我们经常没有足够的时间或数据来正确地做这件事。例如，如果美联储没有在 2008 年后向经济注入这么多资金，我们永远不会知道会发生什么。我们只能观察到发生的一组事件——所以我们不仅错过了红色圆锥(没有美联储印钞的结果)，我们甚至没有观察到绿色圆锥的变化(我们只看到一条绿线)。这就是为什么学者和专家可以争论到面红耳赤，量化宽松和大规模财政刺激是否真正有效，但事实是没有人真正知道(尽管大萧条提供了有用的教训，但它不是 2008 年的可比对照案例)。

事实上，我们只看到一系列事件

结论

如果听起来很累，那是因为确实很累。想得太多太累了！但我们做得越多，它就越成为习惯或根深蒂固的肌肉记忆。它提高了我们的分析能力——最重要的是，它提高了我们独立于过程产生的结果评估过程的能力。

良好的决策(无论是在个人、团体还是公司层面)需要一个强大的、可重复的分析过程。但是每一个决定，无论多么精明，仍然是一个赌注——即使是最有利的赌注也仍然受到偶然事件的影响。通过正确认识机会的作用，我们防止结果扭曲我们的过程(或模型)。通过磨练我们现实地预测采取行动(或不采取行动)的结果范围(以及概率)的能力，我们可以做到即使在不确定的世界中，胜算也是对我们有利的。

方差:样本与总体

原文：https://towardsdatascience.com/variance-sample-vs-population-3ddbd29e498a?source=collection_archive---------26-----------------------

统计 101

ML 中 EDA 的区别以及为什么使用一个而不是另一个

克里斯·贾维斯在 Unsplash 上的照片

在统计学和机器学习中，当我们谈论人口时，我们指的是随机变量的可能值的整个宇宙。如果你知道总体总是可行的，计算平均值和方差如下:

其中 n 是总体的基数(元素的数量)。

大多数时候，我们不能使用整个人口，因为它太复杂了，或者根本不可行。例如，考虑一个问题，当你想分析森林中橡树的高度时。当然，您可以测量森林中的每一棵树，从而收集整个森林的统计数据，但这可能非常昂贵，并且需要很长时间。相反，你可以获取一个样本，比如说，20 棵树，并尝试将样本统计数据和种群统计数据联系起来。因此，对于 N 个样本，我们有:

你现在可以问了:为什么 N-1 而不是 N ？要回答这个问题，我们需要做一些计算。首先，我们计算 s 2 的期望值:

然后，用一点代数知识:

现在，记住:

我们有:

这是什么意思？ 使用 N-1 的公式得出样本方差，平均起来等于未知总体方差。 所以，同样用很少的样本，我们就可以得到一个实际但未知的总体分布参数的合理估计。

如果我们用 N 而不是 N-1 来计算会怎么样？让我们看看:

因此，当我们使用 N 而不是 N-1，时，我们会有一个称为统计偏差的误差，这意味着样本方差(估计量)与真实总体参数(在本例中为方差)存在系统性差异。

校正因子 N/(N-1) 被称为 贝塞尔因子 并允许获得无偏方差 s 的平均值作为有偏方差的函数:

如何在机器学习问题上应用它？当我们试图对一个 ML 问题建模时，我们首先查看数据(在监督和非监督学习中)，搜索模式，统计参数，维度减少，特征选择等。这被称为探索性数据分析或 EDA 。无论是数值特征还是分类特征，我们首先要看观察值的分布。所以，最初，我们必须估计两个参数:平均值和方差。如果 N ，即观察值的数量很小，我们必须应用贝塞尔因子来更好地估计真实方差。如果 N 很大(我们会看到有多大)，我们可以省略这一项，因为 N/(N-1)大约等于 1。

所以，在一个问题上:

• 很少观察到
• 要素的值非常少(因为我们有一个有偏差的数据集)

我们必须应用校正因子。

在任何情况下，我们都不能对结果有信心，因为我们使用的是样本而不是总体。我们所能做的最好的事情是估计一个范围的值，在这个范围内实际方差落在(总体方差的置信区间)内。

让我们看一个例子。想象一片有 10000 棵橡树的森林:这是全部人口。我们想估计高度的分布。假设我们不知道身高呈正态分布，平均值为 10m，标准差(方差的平方根)为 2m。这些是整个人口的统计参数。我们试图通过随机抽取的 20 棵橡树样本来估算这些值。我们重复这个实验 100 次。以下是 Python 中的代码:

import numpy as np
import pandas as pd
from random import sample,choice
import matplotlib.pyplot as plt
import seaborn as sns
import scipynp.random.seed(1234)
mu = 10
sigma = 2
pop = np.random.normal(mu, sigma, 10000)
count, bins, ignored = plt.hist(pop, 100, density=True, color = 'lightgreen')
sns.set_style('whitegrid')
tests = 100
sam = []
mean =[]
std_b = []
std_u = []
fig, axs = plt.subplots(ncols=2)
sns.kdeplot(pop, bw=0.3, ax=axs[0])
for i in range(tests):
    sam_20 = np.random.choice(pop, 20)
    sns.kdeplot(sam_20, bw=0.3, ax=axs[1])
    sam.append(sam_20)
    mean.append(np.mean(sam_20))
    std_b.append(np.std(sam_20))
    std_u.append(np.std(sam_20, ddof=1))frame = {'mean':mean, 'std_b': std_b, 'std_u': std_u}
table = pd.DataFrame(frame)

这样，我们为 100 个实验*中的每一个实验获得样本均值"*均值"标准差(有偏)" std_b" 和无偏( ddof=1 与 (N-1) 相同，而不是 N )标准差" std_u" 。从图形上看，我们得到:

假设 20 个样本是 10000 个总体项目的一个很小的子集，每次我们运行测试，我们都会得到不同的分布。无论如何，平均起来，我们得到了真实均值和标准差的合理估计。我们可以在 Python 中使用以下命令看到这一点:

table.describe()

如您所见，平均而言，无偏样本标准差比有偏样本标准差更接近总体参数的值。

通过置信度分析来衡量这些评估的质量，我们将在下一次讨论。

结论

在本文中，我们已经看到了样本方差受统计偏差的影响，这是由于与整个总体的基数相比，很少观察到的数据发生了失真。我们学习了一种通过贝塞尔因子来减少这种误差的方法，并给出了一个例子。

可变自动编码器

原文：https://towardsdatascience.com/variational-autoencoders-63191b75c576?source=collection_archive---------17-----------------------

VAE 和在哪里找到他们

因此，在过去两年的机器学习面试中，我被要求解释变分自动编码器(VAE)大约五次。对于某个公司来说，这个问题实际上是在同一天由两个不同的面试官重复提出的。

为什么 VAE 是如此受欢迎的面试问题？嗯，它们确实涵盖了概率论和信息论中的一些有趣的主题，而且它们被广泛使用，被认为是相对普通的知识。这些模型有一个独特的规格和一个长期的衍生。这意味着他们可以花 30 分钟走一遍，但仍然很容易让面试官验证。此外，人们通常只了解 VAE 背后的直觉，所以这是一个很好的面试问题(在我看来)。

显然，我真的不喜欢把 VAE 作为一个面试问题，但我也认为他们的申请由于其他原因被夸大了。在以后的文章中，我将推导出我所看到的这类“优化游戏”的真正的潜在机制。现在，让我们彻底了解一下 VAE，作为 VAE 召回的一站式商店。

自动编码器

首先，自动编码器是神经网络的一种形式，其使用一些损失函数 L(x，r(x))专门训练重构函数 r = g(f (x))，以便产生更多信息和/或更低维度的数据重构。这个重构函数是编码函数 z = f(x，θ)和解码函数 x = g(z，θ)的组合，编码函数将数据映射到学习的潜在空间，解码函数 x = g(z，θ)充当生成函数，在潜在空间中为 z 产生 x。这种重建是通过无监督学习创建的，并且通常独立于任何标记或预测任务。

虽然神经网络最常用于生成高维参数化函数以解决判别任务，但自动编码器被训练来重新生成模型被训练的数据。如果我们对我们的损失函数进行正则化，并对我们的网络进行战略性结构调整，我们可以改变自动编码器的行为。

通过引入正则化子ω，我们可以影响我们的自动编码器学习的空间类型。一种这样的增强是稀疏自动编码器(SAE ),它学习稀疏潜在表示。考虑正则化损失函数项 L(x，r(x))+ω(z)，其中ω(z)=λ∑| z _ I |惩罚我们的潜在表示的分量的幅度。这种损失使自动编码器的编码函数 z = f (x，θ)发生偏移，以根据λ的标度为潜在表示的尽可能多的分量产生零或小的值。

作为正则化损失函数的一种替代方法，我们也可以考虑对输入到自动编码器的数据进行正则化。去噪自动编码器(DAE)旨在通过最小化准则 L_DAE = E[L(x，r(N(x))]]来生成简化的重建 r(x)，其中对于某些高斯误差ε = N (0，σ)，N(x) = x + ε。输入到网络中的数据本质上是被破坏的，模型学习过滤掉噪声，以便再现源数据。[1]

这些模型使用什么类型的信息来产生这些潜在的表征？自动编码器学习训练数据的概率密度 p(x)的流形结构
∂对数 p(x)/∂x。这种流形结构有助于引导数据潜在空间的构建。形式上，对于去噪自动编码器，可以示出我们产生以下最优重构函数 r*(x) = x + σ^2 ∂对数 p(x)/∂x + o(σ),其行为渐近于高斯噪声的方差。渐近行为(r *(x)-x)/σ→∂对数 p(x)/∂x 为
σ → 0 允许我们通过 p(x)上的诱导流形确定 x 的信息表示。[2]

可变自动编码器

变分自动编码器(VAE)产生了我们的潜在代表符合给定分布的愿望，以及变分推理过程的简单近似使计算易处理的观察。我们将回顾定义潜在空间分布的步骤，以及使用变分推论来训练我们的模型的步骤。

VAE 的嵌入机制是一种学习的编码分布 z = q(x ),可以更好地对数据集的潜在结构进行排序。根据这种分布，可以构建概括的推理模型，以便重现数据。通过概括，我们意味着当通过给定的生成函数 g(z)馈送时，从在潜在空间 Z 中的样本 Z 上学习的分布 P(Z)中抽取的样本更有可能产生我们的原始数据集 X 的样本。母函数 g(z)可以被视为参数化函数 g(z，θ)或数据 P(X | z；θ).形式上，当我们构建我们的潜在空间时，我们试图最大化观察每个 X 的概率，用于由下式定义的生成过程:

为了评估这个函数，我们需要建立假设的先验分布 P(Z ),以及训练的分布 P(X | Z；θ).先验分布通常假设为正态分布，使得 P(Z) = N(0，I)。

对于 VAE，解码网络用于生成 P(X|z，θ)的样本。我们假设的输出服从(通常为多元正态)分布，均值为 g(z，θ),方差为σ2∫I，如下所示

为了评估这些样本有多精确，我们需要建立一些简化的方法来评估变分推断积分。显然，在 g 的每一个优化步骤中，在整个区域内对样本进行充分的积分计算是非常昂贵的。

为了简化对 P(X|z)的评估，我们只关注(近似地)对那些可能产生 X 的 z 进行采样，服从生成过程 g。我们选择用来表示该采样的分布是 Q(z|X)，由我们的编码函数 f 实现。然后我们用下式近似 P(X|z)

现在，我们需要把期望 E[P(X|z)]和 P(X)联系起来。我们可以通过对潜在空间上任意分布φ(z)的 KL 散度和交叉熵公式进行一些代数操作来做到这一点。

从 KL(φ(z)| | P(z | X))，我们可以导出(在一些对数、互信息和交叉熵替换之后)

如果我们设置φ(z)= Q(z | X)，我们可以看到 log P(X)由我们的编码器 f∞Q(z | X)、我们的解码器 g∞P(X | z)和我们假设的先验分布 P(z)决定。变分自动编码器的主要工程发展是可以有效地优化上述等式的右边。如果 Q(z|X)和 P(z)都假定为高斯分布，KL(Q(z|X)||P(z))有一个封闭形式的解。我们可以引入函数μ(X，θ)和σ(X，θ)，作为编码神经网络层，这又允许我们指定 Q(z|X) = N(μ(X)，σ(X))。VAE 损耗方程的闭合形式的 KL 发散项 KL(Q(z|X)||P(z))是[3]

我们如何强制我们的编码分布 Q(z|X)遵循正态分布 N(μ(X)，σ(X))的约束？在训练期间从
Q(z|X) = μ(X，θ)+εσ(X，θ)采样，噪声项为ε ~N(0，1)，允许我们强制我们的样本表现为 Q 上的期望正态分布。这被称为*重采样技巧。*在评估期间，我们只需要设置ε=0，我们的输出就是确定性的。

结尾部分

现在我们有了 VAE 损失函数的推导，以及 VAE 建筑核心力学的概率解释的基本内容。

在实践中，特别是对于更复杂的层，实际上很难证明由变分体系结构产生的表示实际上遵循其公式中使用的先验分布。这部分是因为在进行任何类型的降维时都会丢失信息。然而，我也认为大多数类似 VAE 的模型在学术数据集之外不太适合他们的先验知识。

在任何情况下，我都会在以后的文章中，试着回顾一下我对变分步骤的更多理论上的抱怨。

参考

[1]伊恩·古德费勒，约舒阿·本吉奥和亚伦·库维尔。*深度学习。*麻省理工学院出版社。2016.

[2]纪尧姆·阿兰和约舒阿·本吉奥。正则化自动编码器从数据生成分布中学习到什么。2014 年

[3]卡尔·多尔什。变型自动编码器教程。 2016 年。

使用 Keras 作为生成模型的可变自动编码器

原文：https://towardsdatascience.com/variational-autoencoders-as-generative-models-with-keras-e0c79415a7eb?source=collection_archive---------8-----------------------

使用 Keras 的可变自动编码器和图像生成|图片由“我的镜头生活”来自 Unsplash | 图片来源

实践教程

引入变分自动编码器和使用 Keras 生成图像

概观

本文旨在让读者对变型自动编码器有一些基本的了解，并解释它们在机器学习和人工智能方面与普通自动编码器的不同之处。不同于普通的自动编码器(如稀疏自动编码器、去噪自动编码器。etc)，变分自动编码器(VAEs)是类似 GANs ( 生成对抗网络)的生成模型。这篇文章主要集中在变化的自动编码器上，我将很快在我即将发布的帖子中写关于生成敌对网络的内容。

在本教程中，我们将讨论如何从头开始用 Keras(TensorFlow，Python)训练一个变分自动编码器(VAE)。我们将以一个简单 VAE 的生成能力的演示来结束我们的研究。

本教程的其余内容可分为以下几类-

1. 背景:变分自动编码器
2. 在喀拉斯建造 VAE
3. 在 MNIST 数据集上训练 VAE
4. 结果
5. 图像生成能力
6. 汇总
7. 延伸阅读和资源

1.变分自动编码器

背景

自动编码器基本上是一种神经网络，它将高维数据点作为输入，将其转换为低维特征向量(即潜在向量)，并且稍后仅利用潜在向量表示来重构原始输入样本，而不会丢失有价值的信息。任何给定的自动编码器都由以下两部分组成——编码器和解码器。模型的编码器部分获取输入数据样本，并将其压缩成潜在向量。而解码器部分负责从学习的(由编码器在训练期间学习的)潜在表示中重建原始输入样本。要了解更多基础知识，请查看我关于 Keras 中的自动编码器和深度学习的文章。

让我们继续考虑，直到现在我们都在同一页上。

普通自动编码器的一个问题是它们独立地编码每个输入样本。这意味着属于同一类的样本(或属于同一分布的样本)可能学习到非常不同的(潜在空间中的远距离编码)潜在嵌入。理想情况下，同类的潜在特征应该有些相似(或者潜在空间上更接近)。这是因为我们没有明确强制神经网络学习输入数据集的分布。由于这个问题，我们的网络可能不太擅长重建相关的看不见的数据样本(或不太可归纳)。

在过去关于 Keras 中的自动编码器和深度学习的教程中，我们训练了一个普通的自动编码器，并学习了 MNIST 手写数字图像的潜在特征。当我们在具有相应标签的潜在空间中绘制这些嵌入时，我们发现相同类别的学习嵌入有时出现相当随机，并且在不同类别的嵌入聚类之间没有清晰可见的边界。下图显示了分布情况-

使用 Keras 的可变自动编码器和图像生成|作者图片| 图片来源

可变自动编码器

变分自动编码器本质上略有不同。它不是直接从输入样本中学习潜在特征，而是实际学习潜在特征的分布。假设输入数据的潜在特征遵循标准正态分布。这意味着学习到的潜在向量应该是以零为中心的，并且它们可以用两个统计量来表示——均值和方差(因为标准正态分布可以仅用这两个统计量来表示)。

因此，变分自动编码器(vae)计算每个样本的潜在向量的均值和方差(而不是直接学习潜在特征),并迫使它们遵循标准的正态分布。因此，网络的瓶颈部分用于学习每个样本的均值和方差，我们将定义两个不同的全连接(FC)层来计算这两个值。VAEs 确保潜在空间中彼此非常接近的点代表非常相似的数据样本(相似的数据类别)。我们将在本教程中证明这一事实。

在进入实现细节之前，让我们首先了解一下 KL-divergence，它将被用作我们模型中的两个优化度量之一。

库尔贝克-莱布勒(KL)散度

在上一节中，我们讨论了对输入数据集的潜在要素实施标准正态分布。这可以使用 KL-散度统计来实现。 KL-divergence 是两个概率分布之间差异的统计度量。因此，我们将利用 KL-散度值作为目标函数(连同重建损失),以确保学习分布非常类似于真实分布，我们已经假设真实分布是标准正态分布。

在这种情况下，最终目标可以写成-

*Objective = Reconstruction Loss + KL-Loss*

这里，重建损失项将鼓励模型学习正确重建原始图像(如果不完全相同，则为同一类别的图像)所需的重要潜在特征。而 KL 散度损失项将确保学习的分布类似于真实分布(标准正态分布)。这进一步意味着分布以零为中心，并且在空间中分布均匀。我们将在教程的后半部分证明这一点。

2.在喀拉斯建造 VAE

最后一节解释了机器学习(ML)和人工智能(AI)中可变自动编码器(VAEs)背后的基本思想。在这一节中，我们将使用 Python 中的 Keras 构建一个卷积变分自动编码器。该网络将在 Keras 数据集中可用的 MNIST 手写数字数据集上进行训练。

为逐步理解和简化起见，本节可分为以下几个部分

1. 数据准备
2. 建筑编码器
3. 潜在分布和采样
4. 建筑解码器
5. 建筑 VAE
6. 损失

数据准备

在本节中，我们将下载 MNIST 手写数字数据集并将其加载到我们的 Python 笔记本中，以开始准备数据。

这是预先加载的依赖项-

**import** numpy as np
**import** matplotlib.pyplot as plt
**import** pandas as pd
**import** seaborn as sns
**import** warningswarnings.filterwarnings('ignore')
**%**matplotlib inline

以下 python 代码可用于下载 MNIST 手写数字数据集。下面还显示了一些示例图像-

**from** tensorflow.keras.datasets **import** mnist(trainX, trainy), (testX, testy) **=** mnist.load_data()print('Training data shapes: X=%s, y=%s' **%** (trainX.shape, trainy.shape))print('Testing data shapes: X=%s, y=%s' **%** (testX.shape, testy.shape))**for** j **in** range(5):
    i **=** np.random.randint(0, 10000)
    plt.subplot(550 **+** 1 **+** j)
    plt.imshow(trainX[i], cmap**=**'gray')
    plt.title(trainy[i])
plt.show()

数据集已经分为定型集和测试集。训练数据集具有分辨率为 28*28 的 60K 手写数字图像。尽管测试数据集由具有相似尺寸的 10K 手写数字图像组成

MNIST 数据集|使用 Keras 的可变自动编码器和图像生成

数据集中的每个图像都是一个 2D 矩阵，表示从 0 到 255 的像素强度。我们将首先归一化像素值(使它们在 0 和 1 之间)，然后为图像通道添加一个额外的维度(由 Keras 的 Conv2D 层支持)。以下是 python 中的预处理代码-

train_data **=** trainX.astype('float32')**/**255
test_data **=** testX.astype('float32')**/**255train_data **=** np.reshape(train_data, (60000, 28, 28, 1))
test_data **=** np.reshape(test_data, (10000, 28, 28, 1))print (train_data.shape, test_data.shape)Out[1]: (60000, 28, 28, 1) (10000, 28, 28, 1)

建筑编码器

在这一节中，我们将定义 VAE 模型的编码器部分。当输入数据类型是图像时，自动编码器的编码器部分通常包括多个重复的卷积层，然后是池层。变型自动编码器的编码器部分也非常相似，只是瓶颈部分与上面讨论的略有不同。

这里是编码器部分的 python 实现，带有 Keras-

**import** tensorflowinput_data **=** tensorflow.keras.layers.Input(shape**=**(28, 28, 1))encoder **=** tensorflow.keras.layers.Conv2D(64, (5,5), activation**=**'relu')(input_data)encoder **=** tensorflow.keras.layers.MaxPooling2D((2,2))(encoder)encoder **=** tensorflow.keras.layers.Conv2D(64, (3,3), activation**=**'relu')(encoder)encoder **=** tensorflow.keras.layers.MaxPooling2D((2,2))(encoder)encoder **=** tensorflow.keras.layers.Conv2D(32, (3,3), activation**=**'relu')(encoder)encoder **=** tensorflow.keras.layers.MaxPooling2D((2,2))(encoder)encoder **=** tensorflow.keras.layers.Flatten()(encoder)encoder **=** tensorflow.keras.layers.Dense(16)(encoder)

上面的片段压缩了图像输入，并将其降低到 16 值特征向量，但这些不是最终的潜在特征。下一节将通过添加潜在特征计算逻辑来完成编码器部分。

潜在分布和抽样

该部分负责从上一部分中提取复杂的特征，并计算潜在特征的平均值和对数方差(因为我们已经假设潜在特征遵循标准正态分布，并且该分布可以用平均值和方差统计值来表示)。两个独立的全连接(FC 层)层用于计算给定数据集的输入样本的平均值和对数方差。

标准正态分布(SND)的这些属性(平均值和对数方差)随后用于估计相应输入数据点的潜在编码。下面定义的函数 sample_latent_features 取这两个统计值并返回一个潜在编码向量。该潜在编码被传递给解码器作为图像重建目的的输入。

def sample_latent_features(distribution):
    distribution_mean, distribution_variance = distribution
    batch_size = tensorflow.shape(distribution_variance)[0]
    random = tensorflow.keras.backend.random_normal(shape=(batch_size, tensorflow.shape(distribution_variance)[1]))
    return distribution_mean + tensorflow.exp(0.5 * distribution_variance) * random

distribution_mean = tensorflow.keras.layers.Dense(2, name='mean')(encoder)distribution_variance = tensorflow.keras.layers.Dense(2, name='log_variance')(encoder)latent_encoding = tensorflow.keras.layers.Lambda(sample_latent_features)([distribution_mean, distribution_variance])

这些潜在特征(从学习的分布计算)实际上完成了模型的编码器部分。现在编码器模型可以定义如下

encoder_model **=** tensorflow.keras.Model(input_data, latent_encoding)encoder_model.summary()

编码器非常简单，只有大约 57K 个可训练参数。

建筑解码器

该模型的编码器部分将图像作为输入，并将该图像的潜在编码向量作为输出，该输出是从输入数据集的学习分布中采样的。解码器的工作是将该嵌入向量作为输入，并重新创建原始图像(或与原始图像属于相似类别的图像)。由于潜在向量是特征的非常压缩的表示，解码器部分由多对去卷积层和上采样层组成。去卷积层基本上与卷积层相反。上采样层用于恢复图像的原始分辨率。通过这种方式，它可以重建原始尺寸的图像。

下面是用 TensorFlow 的 Keras API 实现的解码器部分的 python 实现

decoder_input = tensorflow.keras.layers.Input(shape=(2))
decoder = tensorflow.keras.layers.Dense(64)(decoder_input)
decoder = tensorflow.keras.layers.Reshape((1, 1, 64))(decoder)
decoder = tensorflow.keras.layers.Conv2DTranspose(64, (3,3), activation='relu')(decoder)

decoder = tensorflow.keras.layers.Conv2DTranspose(64, (3,3), activation='relu')(decoder)
decoder = tensorflow.keras.layers.UpSampling2D((2,2))(decoder)

decoder = tensorflow.keras.layers.Conv2DTranspose(64, (3,3), activation='relu')(decoder)
decoder = tensorflow.keras.layers.UpSampling2D((2,2))(decoder)

decoder_output = tensorflow.keras.layers.Conv2DTranspose(1, (5,5), activation='relu')(decoder)

解码器模型对象可以定义如下

decoder_model **=** tensorflow.keras.Model(decoder_input, decoder_output)decoder_model.summary()

解码器同样简单，具有 112K 可训练参数。

构建 VAE(可变自动编码器)

最后，可以通过组合编码器和解码器部分来定义变分自动编码器(VAE)。下面是你如何通过在编码器后面粘贴解码器来创建 VAE 模型对象。

encoded **=** encoder_model(input_data)decoded **=** decoder_model(encoded)autoencoder **=** tensorflow.keras.models.Model(input_data, decoded)autoencoder.summary()

整个设置非常简单，只有 170，000 个可训练模型参数。写目标(或优化函数)函数的时间。

失败

如前所述，变分自动编码器(VAE)的最终目标(或损失)函数是数据重建损失和 KL 损失的组合。在本节中，我们将通过结合这两个统计数据来定义我们的自定义损失。

get_loss 函数的以下实现返回一个 total_loss 函数，它是重建损失和 KL-loss 的组合，定义如下

Github 链接

最后，让我们编译模型，为培训做好准备-

autoencoder.compile(loss=get_loss(distribution_mean, distribution_variance), optimizer='adam')

3.训练 VAE(可变自动编码器)

就像普通的自动编码器一样，我们将通过输入和输出完全相同的图像来训练它。该模型被训练 20 个时期，批次大小为 64。

这是训练总结-

autoencoder.fit(train_data, train_data, epochs=20, batch_size=64, validation_data=(test_data, test_data))

我希望它可以被训练得更好一点，但这是验证损失没有太大变化的地方，我继续进行。

4.结果

在本节中，我们将在测试图像上看到我们的模型的重建能力。以下 python 脚本将从测试数据集中选取 9 幅图像，我们将为它们绘制相应的重建图像。

offset**=**400
print ("Real Test Images")# Real Images**for** i **in** range(9):
    plt.subplot(330 **+** 1 **+** i)
    plt.imshow(test_data[i**+**offset,:,:, **-**1], cmap**=**'gray')
plt.show()# Reconstructed Imagesprint ("Reconstructed Images with Variational Autoencoder")**for** i **in** range(9):
    plt.subplot(330 **+** 1 **+** i)
    output **=** autoencoder.predict(np.array([test_data[i**+**offset]]))
    op_image **=** np.reshape(output[0]*****255, (28, 28))
    plt.imshow(op_image, cmap**=**'gray')
plt.show()

这是输出-

结果|使用 Keras 的可变自动编码器和图像生成

上述结果证实了该模型能够以相当高的效率重建数字图像。然而，这里要注意的一件重要事情是，一些重建图像在外观上与原始图像非常不同，而类别(或数字)总是相同的。发生这种情况是因为重建不仅仅依赖于输入图像，它是已经学习的分布。并且这种学习到的分布是模型输出中引入变化的原因。这很有趣，不是吗！

这里要注意的第二件事是输出图像有点模糊。这是变化的自动编码器的常见情况，由于潜在向量(瓶颈)非常小，并且如前所述，存在学习潜在特征的单独过程，它们经常产生有噪声的(或质量差的)输出。变分自动编码器(vae)实际上不是为了重建图像而设计的，真正的目的是学习分布(这给了他们生成假数据的超能力，我们将在后面的帖子中看到)。

潜在特征聚类

如前所述，变分自动编码器(vae)学习潜在特征的潜在分布，这基本上意味着属于同一类的样本的潜在编码在潜在空间中不应彼此相距太远。其次，总体分布应该是标准正态的，应该是以零为中心。

让我们为所有的测试图像生成潜在嵌入，并绘制它们(相同的颜色代表属于相同类别的数字，取自地面真实标签)。这是 python 代码-

x = []
y = []
z = []for i in range(10000):
    z.append(testy[i])
    op = encoder_model.predict(np.array([test_data[i]]))
    x.append(op[0][0])
    y.append(op[0][1])df = pd.DataFrame()
df['x'] = x
df['y'] = y
df['z'] = ["digit-"+str(k) for k in z]plt.figure(figsize=(8, 6))
sns.scatterplot(x='x', y='y', hue='z', data=df)
plt.show()

基于 Keras 的变分自动编码器和图像生成

上图显示分布以零为中心。相同类别数字的嵌入在潜在空间中更接近。数字分离边界也可以很容易地画出来。这正是我们想要从可变自动编码器中实现的。让我们跳到最后一部分，测试我们模型的生成能力。

5.伪图像生成

可变自动编码器可以用作生成模型。前面的部分显示了输入数据的潜在编码遵循标准的正态分布，并且对于不同类别的数字存在明显的界限。

请想一想——如果我们已经知道，空间的哪一部分用于哪个类，我们甚至不需要输入图像来重建图像。这意味着我们实际上可以通过从空间(潜在分布空间)传递随机点来生成与训练数据集具有相似特征的数字图像。以这种方式，变分自动编码器可以用作生成模型，以便生成假数据。

正如我们所看到的，潜在编码的范围在[-3 到 3 之间，在 x 轴上也是-3 到 3]。让我们用只属于这个范围的随机潜在编码生成一串数字。

generator_model **=** decoder_modelx_values **=** np.linspace(**-**3, 3, 30)
y_values **=** np.linspace(**-**3, 3, 30)figure **=** np.zeros((28 ***** 30, 28 ***** 30))**for** ix, x **in** enumerate(x_values):
    **for** iy, y **in** enumerate(y_values):
        latent_point **=** np.array([[x, y]])
        generated_image **=** generator_model.predict(latent_point)[0]
        figure[ix*****28:(ix**+**1)*****28, iy*****28:(iy**+**1)*****28,] **=** generated_image[:,:,**-**1]plt.figure(figsize**=**(15, 15))
plt.imshow(figure, cmap**=**'gray', extent**=**[3,**-**3,3,**-**3])
plt.show()

生成的数据|使用 Keras 的可变自动编码器和图像生成

你可以在上面的图像矩阵中找到所有的数字(从 0 到 9 ),因为我们试图从潜在空间的所有部分生成图像。生成变化的笔迹的能力是不是太棒了！

6.摘要

本教程解释了深度学习和人工智能中的可变自动编码器。通过一个基本的介绍，展示了如何用 python 中的 Keras 和 TensorFlow 实现 VAE。在 MNIST 手写数字数据集上进一步训练该模型，并给出重建结果。

我们已经看到，潜在的编码遵循标准的正态分布(全部归功于 KL-divergence ),以及如何将模型的训练解码器部分用作生成模型。我们已经通过仅使用模型的解码器部分生成假数字来证明了这一说法。

如果你有兴趣阅读我关于去噪自动编码器的文章

卷积去噪图像降噪自动编码器

Github 代码链接:https://github.com/kartikgill/Autoencoders

原载于 滴艾 。

感谢阅读！希望这对你有帮助。请在下面评论，让我知道你的反馈。下一篇文章再见。

阅读下一篇>>>

1. Keras 和深度学习中的自动编码器(简介)
2. 优化者解释训练神经网络
3. 用量化技术优化张量流模型
4. 用 PyTorch 进行深度学习:简介
5. 用 PyTorch 进行深度学习:第一个神经网络

参考

1. 如何在 Keras 中构建一个可变的自动编码器
2. 变型自动编码器
3. 【https://keras.io/examples/generative/vae/

相关研究论文

1. 语法变分自动编码器
2. 用于分子图生成的连接树变分自动编码器
3. 图像、标签、字幕深度学习的变分自动编码器
4. 使用重构概率的基于变分自动编码器的异常检测
5. 用于文本生成的混合卷积变分自动编码器

用于假人的可变自动编码器(VAEs)——循序渐进教程

原文：https://towardsdatascience.com/variational-autoencoders-vaes-for-dummies-step-by-step-tutorial-69e6d1c9d8e9?source=collection_archive---------3-----------------------

DIY 实践指南与实践代码的建设和培训与 Keras 的名人脸上的 VAEs

本文介绍了您需要从生成模型中获得的一切。我们提供了如何在大型图像数据集上训练条件值并使用它们生成新的标记图像的分步指南。

Pladicon 在 Pixabay 拍摄的照片

动机

既然已经有这么多数据，为什么我们还需要生成新的数据？据 IDC 称，目前有超过 18 兆字节的数据。

大多数机器学习任务都需要标记数据。获取高质量、有标签的数据是困难的。如果我们自己生成这些数据，我们就可以随心所欲地利用它们。新的数据可以给我们想法和选择。

如何生成看不见的图像？

在这篇文章中，你将学习什么是变分自动编码器，以及如何创建自己的变分自动编码器来生成新的和看不见的图像。我们不用数学来解释潜在的概念和直觉。

使用我们的 VAE 码(自创)的图像及其重建的例子

数据

我们将使用知名名人数据集的子集来帮助我们建立面部生成模型。数据集可以按照名人脸网站上的描述下载。它提供了一个大规模的人脸属性数据集，包含超过 20 万张名人图片，每张图片都有 40 个属性注释。

• 10177 个身份，
• 202，599 个面部图像，
• 5 个地标位置，以及
• 每个图像 40 个二进制属性注释。

下面我们从名人数据集中随机挑选一些面孔，并显示他们的元数据(属性)。图像高度为 218 像素，宽度为 178 像素，有 3 个颜色通道。

什么是自动编码器(AE)？

通过观察数以千计的名人面孔，神经网络可以学习生成不存在的人的面孔。

神经网络是我们可以用来获得函数近似值的许多可能方法之一。他们受欢迎的一个原因是他们学习表达的能力。假设我们为网络提供了正确的标签，网络可以学习在将图像分类为狗或猫时重要的特定表示。这是监督学习。

在某些情况下，我们没有这些标签。然而，我们可以训练两个网络，一个学习表示，另一个通过最小化重建损失函数从表示中重建。这是一个自动编码器。它之所以有这个名字，是因为它会自动找到编码输入的最佳方式，以便解码后的版本尽可能接近输入。

自动编码器由一对两个相连的神经网络组成:编码器模型和解码器模型。它的目标是找到一种方法将名人的脸编码成一种压缩的形式(潜在空间)，以这种方式重建的版本尽可能接近输入。

自动编码器的工作组件(自行创建)

编码器模型将输入 x 转换为小型密集表示 z ，类似于卷积神经网络通过使用滤波器来学习表示的工作方式。

解码器模型可以被视为能够生成特定特征*x’*的生成模型。

编码器和解码器通常作为一个整体来训练。损失函数惩罚网络创建不同于输入面的输出面。

因此，编码器学会在有限的潜在空间中保存尽可能多的相关信息，并巧妙地丢弃不相关的部分，例如噪声。

解码器学习获取压缩的潜在信息，并将其重建为完整的名人脸。

自动编码器也可以用于维度减少和图像去噪，但也可以成功用于无监督机器翻译。

什么是变分自动编码器(VAE)？

通常，编码器产生的潜在空间 z 人口稀少，这意味着很难预测该空间中值的分布。在 2D 制图表达中，值是分散的，空间似乎得到了很好的利用。

这对于压缩系统来说是一个非常好的特性。然而，对于生成新的名人图像，这种稀疏性是一个问题，因为找到解码器将知道如何产生有效图像的潜在值几乎是不可能的。

此外，如果空间在簇之间有间隙，并且解码器从那里接收到变化，它将缺乏生成有用的东西的知识。

变分自动编码器通过使潜在空间更可预测、更连续、更少稀疏来工作。通过强制潜在变量成为正态分布，VAEs 获得对潜在空间的控制。

从 AE 到 VAE 使用随机变量(自创)

VAEs 不是将潜在值直接转发给解码器，而是使用它们来计算平均值和标准偏差。然后从相应的正态分布中对解码器的输入进行采样。

在训练过程中，VAEs 通过在损失函数中包含 kull back-lei bler 散度，迫使该正态分布尽可能接近标准正态分布。VAE 将会改变或探索脸部的变化，而且不只是以随机的方式，而是朝着一个期望的、特定的方向。

条件变分自动编码器允许基于潜在变量 z 和附加信息(如面部元数据(微笑、眼镜、肤色等))对输入进行建模。).

图像数据生成器

让我们建立一个(有条件的)可以在名人脸上学习的 VAE。我们使用定制的 Keras 内存高效生成器来处理我们的大型数据集(202599 张图像，ca。每个 10KB)。这背后的想法是在训练过程中即时获得批量图像。

VAE 网

我们希望编码器是一个卷积神经网络，它获取图像并输出分布参数 Q(z|[x，c]) 其中 x 是人脸的输入图像， c 是条件变量(人脸的属性)，而 z 是潜在变量。对于本文，我们使用一个由两个卷积层和一个池层组成的简单架构。

解码器是一个反过来构建的卷积神经网络。它是一个生成网络，向似然分布 P([x，c]|z) 输出参数。

整个 VAE 网络的架构创建如下。

培养

下图显示了 celebA 数据集中图像的 VAE 模型的学习过程。代码在使用 1 个 GPU 的 AWS 实例上运行了大约 8 个小时。

可视化潜在表征

经过训练后，我们现在可以从我们的数据集中随机选取一幅图像，并使用经过训练的编码器来创建图像的潜在表示。

使用 16 个实数的向量的潜在表示，我们可以可视化解码器如何重建原始图像。

虽然重建的图像很模糊，但我们可以看到与原始图像有很强的相似性:性别、衣服颜色、头发、微笑、肤色。

生成看不见的面孔

条件值能够改变潜在空间以产生新数据。具体来说，我们可以使用解码器生成一个不可见图像的随机计数，同时根据给定的属性对其进行调节。

虽然我们的变分自动编码器产生了模糊和非真实感的人脸，但我们可以识别那些从未存在过的人类的性别、肤色、微笑、眼镜、头发颜色。

给人们一个微笑

条件值可以在属性之间进行插值，使一张脸微笑，或者在以前没有眼镜的地方添加眼镜。下面，我们从数据集中随机选择一张名人脸，并利用潜在表征的改变，将其从女性脸变形为男性脸。我们也改变这些脸来展示原本不存在的微笑。

结论

在本文中，我们介绍了条件变分自动编码器，并演示了它们如何学习如何生成新的标记数据。

我们提供了 Python 代码，用于在大型名人图像数据集上训练 VAEs。该方法和代码可以扩展到多个其他用例。

生成敌对网络(GANs)往往会产生更好看的图像，因为它们学会了区分什么对人类来说是真实感的，什么不是。

[## 虚拟生成对抗网络(GAN)——循序渐进教程

用防弹 Python 代码理解、构建和训练 GANs 的终极初学者指南。

towardsdatascience.com](/generative-adversarial-network-gan-for-dummies-a-step-by-step-tutorial-fdefff170391)

应认真考虑使用和 GAN 技术制作虚假图像、视频或新闻的道德问题，应负责任地使用这些技术。

伊万诺夫·古德在 Pixabay 上的照片

非常感谢 Vincent Casser 的精彩代码，在他的博客中提供了一种更先进的实现图像处理卷积自动编码器的方法。我获得了文森特的明确授权，为这篇文章改编他的 VAE 代码。从头开始构建一个可用的 VAE 相当棘手。功劳归于文森特·卡瑟。

我推荐阅读 Joseph Rocca 的后续文章:了解变分自动编码器(VAEs) 。

感谢安妮·邦纳的编辑笔记。

保持健康和安全！

变分高斯过程(VGP)——当事物不是高斯分布时该怎么办

原文：https://towardsdatascience.com/variational-gaussian-process-what-to-do-when-things-are-not-gaussian-41197039f3d4?source=collection_archive---------7-----------------------

来自 Pixabay 的图片

人们可能会有这样的印象，贝叶斯方法只有在一切都是高斯的情况下才有效。不是这样的。现代贝叶斯机器学习可以轻松处理非高斯分布。

我将描述高斯过程二元分类模型。它使用高斯先验和伯努利似然。我将解释为什么伯努利分布阻止我们使用贝叶斯规则计算后验概率。我们称这样的后验为难处理的。不管棘手与否，我们都需要后验来做预测。然后，我将提出解决方案——变分推论——来逼近难以处理的后验概率。

二元分类问题

我们有一些训练数据 (X，Y) :

X 为实数值； Y 是作为标签的二进制值。任务是设计一个模型，从输入 X 预测标签 Y 。这是一个二元分类问题。让我们绘制数据，并考虑如何使用高斯过程来构建二元分类器:

在上图中，红叉表示训练数据点。我们意识到，如果我们使用高斯过程回归来对数据建模，就会出现问题:高斯过程回归对连续函数 f(x) 进行建模。我们不能要求这样的函数只返回 1 或 0。高斯过程回归模型可能在某些 x 位置返回 1.1，在其他位置返回 0.3。

这个问题很容易解决。我们使用伯努利分布来模拟二进制输出。伯努利分布只有一个参数，即介于 0 和 1 之间的概率，它描述了事件发生的可能性。为了在这个概率和我们的高斯过程函数 f(x) 之间建立联系，我们定义了函数 g(f(x)) ，它将实值高斯过程函数 f(x) 压缩为范围[0，1]内的一个值。

g 这种方式定义的是 x 的函数。就像当你有 f(x) = x+1 并定义 g(f(x)) = f(x) 时，那么 g 就是 x : g(x) = (x+1) 的函数。

许多挤压函数(检查 sigmoid 函数)都可以工作，我们选择 logit 函数是因为它被广泛使用并且它的数学属性很容易解释:

下图说明了这一想法:

该图有三个部分:

1. 顶部是黑色的高斯过程函数 f(x) 。它在对应标签为 1 的 x 位置经历一些大的实数值(相对于范围[0，1])。并且它在 x 位置经过一些小值，在这些位置对应的标签是 0。
2. 中间部分用蓝色显示了 g(x) 函数。它将 f(x) 函数的值压缩到范围[0，1]内。
3. 底部用红色显示我们的训练数据。

让我们浏览一个示例数据点—图中突出显示的数据点*(y₄=1 x₄=9)*。对于这一点，顶部的读数为 f(9)=3 。在中间部分，挤压功能 g(9) 为:

伯努利分布以 0.98 为参数，表示 Y₄=1 的概率为 0.98 。根据底部，标签 Y₄ 确实是 1 *。*这个模式不错。

顺便说一下，这种使用挤压函数的想法与我们如何构建神经网络分类器非常相似:一个 softmax 函数将最后一个隐藏层中神经元的输出值(全实线中的值)转换为 0 到 1 之间的概率。

受到这个例子的鼓舞，我们正式定义了这个贝叶斯模型。但首先，一些符号。

一些符号

中支持文本中的 Unicode。这让我可以写很多数学下标符号，比如 X₁和 Xₙ.但是我不能写下其他的下标。例如:

所以在文中我会用一个下划线“”来引出这样的下标，比如 *X** 和 *X_1 。

如果一些数学符号在你的手机上显示为问号，请尝试从电脑上阅读这篇文章。这是某些 Unicode 呈现的已知问题。

二元分类模型

贝叶斯模型总是由两部分组成:先验和似然。

院长

顾名思义，高斯过程先验总是多元高斯分布。它模拟了上图顶部的函数。为了刷新你的记忆，看一看理解高斯过程，苏格拉底方式。

GP 先验是两个随机变量向量 f(X) 和 f(X_)* 之间的联合多元高斯分布:

• f(X) 为长度为 n 的随机变量向量， n 为训练数据点的个数。 f(X) 代表训练位置 X 的函数值。
• f(X_)* 是另一个随机变量向量。它表示测试位置 X_，的函数值，即您希望模型进行预测的位置。 f(X_) 的长度与测试位置 X_* 的数量相同。

我用 f 作为 f(X) 的简写。该向量中第个随机变量的 fᵢ ，对应于第个第个训练数据点。同样， f_* 是 f(X_)* 的简写。

GP prior 使用以下核函数 k(x，X’)来定义位置 x 处的两个独立随机变量 f(x) 与位置*X’处的f(X’)之间的相关性(位置集合 x 和X’都可以位于位置集合 X 或 X_)

其中 exp 是指数函数。 l 和 σ 为内核参数。 l 称为长度刻度， σ 称为信号方差。 l 和 σ 是标量。由于它们是我们引入模型的参数，我们也称它们为模型参数。内核函数 k 提到了我们在协方差块 k(X，X) 中的训练位置 X 和在协方差块 k(X_，X_) 中的测试位置 X_* 以及在 k(X_，X)* 中的 X 和 X_* 。 k 不提训练数据 Y 。

简而言之，GP 先验在训练位置引入随机变量 f ，在测试位置引入随机变量 f_* ，并使用核函数 k 来定义任意两个独立随机变量之间的相关性。这种相关性使模型能够通过使用训练位置的随机变量的信息来推断测试位置的随机变量的值，换句话说，进行预测。

尽管 GP prior 是 f 和 f_，*的联合分布，但只有 f 与我们的训练数据有关联。我们可以在导出边际分布 p(f) 之前将多变量高斯边际化规则应用于 GP，只是对于 f. 和几个段落之后，我们将通过似然性建立从训练数据 Y 到 f 的联系。

下面是被边缘化的 f ，其中 K 是 k(X，X)的简称:

或者，以概率密度函数的形式:

与符号 p(f) 相比，符号p(f；l，σ ) 强调这个分布是在随机变量 f 之上；并提及模型参数 l 和 *σ。*它们来自于 K 中提到的指数平方核函数。 det 是矩阵行列式运算符。

注意在本文中，我们不把 l 和 σ 当作随机变量。我们将它们视为模型的标量参数。作为标量，每个模型参数可以取一个单个具体值，如 l=10.4 。我们现在不知道这个值。我将介绍一个参数学习来为 l (以及其他模型参数)找到单个值。这个具体值称为点估计。

还有其他的参数学习方法存在，比如我们可以把 l 和 σ 看成随机变量，可以取一个值的分布。这被称为完全贝叶斯治疗，是未来文章的主题。

可能性

让我们引入一个新的随机变量向量 y(X) ，我将其简称为*y。*我们将我们的观察值 Y 建模为随机变量中的样本 y — 每个 Yᵢ 都是来自相应随机变量 yᵢ.的样本 y 的长度为 n ，其中 n 为训练数据点数。

由于 y 是我们有观测值的随机变量，我们称之为观测随机变量。相反，随机变量 f 和 f_* 没有关联观测值，我们称之为潜在随机变量。

表示为 p(y|f，f_)* 的可能性是一个条件分布。它描述了观察随机变量 y 的概率，给出了 GP 先验的潜在变量 f 和 f_* 。

但是我们不应该设计我们的模型，使得有观测值的随机变量依赖于没有观测值的随机变量。所以我们强制要求 y 不依赖于 f_* 。这通过在以下条件下丢弃不相关的 f_* 来产生可能性 p(y|f) :

由于我们的观察值 Yᵢ 是一个二进制值，我们将其建模为来自伯努利随机变量 yᵢ.的样本yᵢ的概率密度函数为:

其中 g ᵢ 是伯努利分布的唯一参数。 gᵢ 是函数应用 g(fᵢ) 的简称。通过为 yᵢ 插入值 Yᵢ ，我们可以评估 yᵢ 获得该值的概率。所以 y 是伯努利随机变量的向量。

让我们进一步破译伯努利概率密度: gᵢ ，是 *g(fᵢ的简写，*代表观察随机变量 yᵢ 的值为 1 的概率。由于 yᵢ 只能取 0 或 1，那么 yᵢ 取值为 0 的概率就是 1-gᵢ 。上面的表达式是一个简洁的形式，它将这两种情况结合成一个公式。

上述伯努利密度函数建立了从观测值 Yᵢ 到观测随机变量 yᵢ.的联系通过调用函数 g 并以 fᵢ 为参数，我们建立了从 Yᵢ到fᵢ的连接:

该函数具有以下形状，x 轴绘图 fᵢ 和 y 轴绘图g(fᵢ:

让我们来看看这个可能性函数在一个地点的全部细节 Xᵢ:

第(1)行是在给定随机变量 fᵢ 的情况下，在单个位置xᵢ的可能性函数的完整公式，该公式来自我们的 GP prior 。

第(2)行将实际观察值 Yᵢ 插入随机变量 yᵢ.

这个可能提到了我们的训练数据 Yᵢ ，但是没有提到训练数据 Xᵢ 。与之前的 GP 相比，它提到了数据 Xᵢ 而不是 Yᵢ 。而随机变量 fᵢ 是先验和似然之间的桥梁，将 Xᵢ 和 Yᵢ 连接在一起。

我们希望我们的模型给可能性一个高值。所以我们来分析一下，当可能性评估为高值时。让我们来关注一下 Yᵢ=1 时的情况(Yᵢ=0 的情况也类似)。当 Yᵢ=1 时，可能性简化为:

为了使这个数量尽可能高(最大值是 1)，我们需要 fᵢ 有一个更大的平均值，并且希望方差小。这样， fᵢ的期望值就是一个大值。与此同时，上述公式的期望值将更接近于 1，即其最大值。

注意，我们并不关心 fᵢ.的确切平均值只要可能性评估为高概率数字，我们就高兴，不管平均值是 4 还是 7。我们可以想象 fᵢ回归到一些潜在的数据空间。

在我们理解了单个数据点 (Xᵢ，Yᵢ) 的可能性之后，就很容易理解所有训练数据点的可能性了。我们假设观察值相互独立，因此所有观察随机变量 p(y|f) 的可能性是所有 n 个个体可能性的乘积，其中 n 是训练数据点的数量:

这里，似然 p(y|f) 没有提到任何模型参数， l 或 σ 。但一般来说，似然函数可以引入新的模型参数。

现在你已经学会了如何将非高斯分布纳入你的可能性。同样，您可以引入其他非高斯分布，如 student-t 来处理异常值。

后面的

在我们描述了先验 p(f) 和似然 p(y|f) 之后，就到了后验 p(f|y) 的时候了。贝叶斯规则根据先验和似然性给出了后验的概率密度函数:

后面提到了我们的两个模型参数 l，σ 。它们来自 prior p(f) 内部的内核函数 k 。有时我们用 p(f | y；l，σ ) 来强调这些模型参数。

分母叫做边际可能性。这是一个 n 维积分，每维是随机变量向量 f 中的一个随机变量。我们必须首先计算这个积分的结果，然后我们可以评估后验 p(f|y) 作为一个函数，例如，进行预测。

计算积分意味着计算出积分结果的解析表达式。例如，我们算出 ∫ x dx 的结果为 x 。表达式 x 是解析的——我们可以插入一个值 x 来得到积分的结果。

然而，我们不知道如何计算这个积分。积分符号内的函数是 p(y|f) 中的伯努利概率密度函数和 p(f) 中的一个多元高斯概率密度函数的乘积。你不能从你的微积分书中找到一个规则来计算这个积分。人们称之为棘手的整合。

你可能想知道，后验不应该是来自 GP 先验的同一组随机变量 f_* 和 f 的分布吗？所以不是 p(f|y) ，后验应该是 p(f_，f|y)* ，不是？你是对的。但是我们所需要的就是和 p(f|y) 一起工作，因为:

换句话说，我们决定将联合后验概率 p(f_，f|y)* 分解为条件概率 p(f_|f)* 和边际概率 p(f|y) 。通过将多元高斯条件规则应用于 GP 先验，我们知道了 p(f_|f)* 的公式。所以我们只需要用 p(f|y) 就可以了。

变分推理

我们不知道如何推导后验 p(f|y)，的解析表达式，换句话说，我们不知道后验的结构。下一个最好的方法是用另一个分布来近似后验概率，该分布的结构在我们的控制中。所以通过设计，我们知道这个新分布的概率密度函数的解析表达式。

请注意，我们希望使用一个新的分布来逼近后验概率 p(f|y) ，而不是逼近边际可能性 p(y) = ∫ p(y|f) p(f) df。后验概率是对 f 的概率分布。我们需要这个分布来做预测。边际可能性是在贝叶斯规则的右侧引入的一个数，以确保左侧是一个积分为 1 的适当概率密度函数。当我们使用一个新的分布来逼近后验概率 p(f|y) 时，我们自然会考虑边际可能性，因为这个新的分布在设计上是合适的*。*

下面是我对“用另一种分布逼近后验概率”这句话的理解:

1. 这个新分布的概率密度函数覆盖与后验概率相同的一组随机变量 f 。所以这个新的分布可以被插入到后验概率出现的地方。
2. 这种新分布的行为类似于后验分布。这意味着对于一个 f ，新分布返回的概率与真实后验概率返回的概率相似。这样，这个新分布的概率密度函数的形状类似于后验分布。

为了理解短语“新发行版的结构在我们的控制之中”，让我们首先进一步发展这个近似概念。

变分分布 q

我们称我们用来近似后验的分布为变分分布。* 我们用 q(f)到表示它的概率密度函数。 q(f) 应该是什么样子？*

我们希望 q(f) 看起来像后验 p(f|y) ，但是我们不知道后验的结构。我们提出 q 是一个均值 μ 协方差矩阵σ的多元高斯分布如何？**

均值和协方差矩阵完全决定了多元高斯分布的形状。随着它们的改变，多元高斯分布的形状也随之改变。因此我们可以找到 μ 和σ的具体值，使得 q(f) 的形状类似于 p(f|y)的形状。我们使用符号q(f；μ，σ)来强调这两个参数。我们称它们为变分参数。

下图说明了制作形状 q(f；μ，σ)接近于p(f | y；l，σ )* ，其中 f 是单个随机变量，而不是随机变量向量:*

紫色分布呈不规则形状，说明后部的结构我们不知道。绿色分布是钟形的，说明了变分分布，这是一个多元高斯分布。我们调整 μ 和σ的值，使绿色尽可能与紫色重叠，使q(f；μ，σ)近似 p(f|y) 。

我们这里有两个非常重要的假设:

1. 我们假设存在一些具体的变分参数值，使得我们的多元高斯分布q(f；μ，σ)类似于非高斯后验p(f | y；l，σ)。**
2. 我们假设可以通过参数学习为 μ 和σ找到这些具体值中的至少一个。**

这两种假设在现实中都可能是错误的。所以对于某些数据集，我们可能无法找到好的 q(f) 的后验近似。这意味着我们可能找不到一个好的贝叶斯模型来很好地拟合训练数据——机器学习任务失败了。这很糟糕，但它发生了。

顺便说一下，变分推理在现代机器学习中非常重要，是一个你无法逃脱的话题。参见其在时间序列预测和人脸生成中的应用。

我们通过 μ和σ来控制 q 的形状

*在多元高斯建议下，*的概率密度函数 q(f；μ，σ)是:

其中 n 是训练数据点的数量， det 是矩阵行列式运算符。σ⁻是协方差矩阵σ的逆。

这个公式也显示了解析公式的重要性——它可以在插入未知数的具体值后进行评估: f ，这是一个浮点向量， μ ，这是一个浮点向量，σ，这是一个大小为 n×n 的浮点矩阵。当我们有了这三个的一些具体值后，上面的公式计算出一个单一的概率数——多元高斯分布下 f 的概率。**

变分分布 q 为多元高斯的提议现在看起来是纯粹的猜测，但稍后你会看到这种结构有很大的优势。

变分分布q(f；μ，σ)具有固定的结构——全多元高斯分布的结构。“满”意味着σ是没有零的密集矩阵。让我们稍后在参数学习部分回到这个问题。现在对我们来说重要的是，我们可以通过两个变化参数 μ 和σ来改变它的形状。理论上，通过改变 μ 和σ*的值，我们可以得到随机变量向量 f 上任何可能的高斯分布。这就是为什么 q 的结构在我们掌控之中。*

一旦我们决定了 μ 和σ的具体值，我们就有了完全指定的 q ，我们可以用它作为真实后验概率的近似值。但是我们如何决定将哪些具体值赋给 μ 和σ*？显而易见的答案是选择 μ 和σ的值，使得*q(f；μ，σ)尽可能靠近后 p(f|y，l，σ ) 。

我们意识到我们需要一个从分布 q(f) 到后验 p(f|y) 的“接近度”度量。我们需要一种方法来最小化这种亲密度。

kull back-lei bler 散度和证据下界

因为我们希望我们的变分分布q(f；μ，σ)向后方靠近p(f | y；l，σ )* ，我们需要一个他们之间的贴近度。我们可以使用 kull back-lei bler 散度(KL 散度)来定义一个分布与另一个分布的接近程度。变分分布之间的 KL 散度公式q(f；μ，σ)和后验p(f | y；l，σ ) 是:*

在上面的推导中，我将模型参数(核参数和变分参数)涂成红色。

线(1)是 KL 发散的符号。KL 运算符内的双竖线||分隔我们要计算距离的两个分布。

第(2)行定义了 KL 背离的含义。它是两个分布之间的预期对数概率比。如果对于某些 f 、 q(f) 和 p(f|y) 评估为相同的概率，那么它们的比值为 1， log(1) 为 0。所以这个特殊的 f 对总的散度度量没有贡献。另一方面，如果该比率不是 1，则该 f 有助于最终的散度测量。因此 KL 散度是在对数标度中 f 的概率密度之间的平均差。

第(3)行用它的数学计算——积分——代替了期望值。它显示被积函数(插入训练数据 X 和 Y )有 4 个未知数: l，σ，μ，σ，这是全套模型参数。随机变量 f 被积分出来。

我们要的是 q 的形状(f；μ，σ)尽可能接近p(f | y；l，σ ) 尽可能。KL 散度 KL(q(f)||p(f|y)) 编码这种接近度，并且它是模型参数的函数。使 q(f) 尽可能接近 p(f|y) 的方法是最小化关于模型参数的 KL 发散:

你内心一定有个疑问在燃烧——我们不知道后验p(f | y)的解析形式，如何才能最小化提到它的 KL 散度？

为了回答这个问题，我们需要继续操作公式，但是现在为了节省空间，我忽略了推导中的模型参数。请记住，无论我们如何操纵 q(f) 和 p(f|y) 之间的 KL 散度，它始终是四个模型参数的函数。

在第(6)行，我们不写 log p(y) 周围的期望运算，因为 log p(y) 没有提到随机变量 f ，忽略的期望是相对于 f 的。所以期望评估为 log p(y) 。我们对 log p(y) 有一个著名的名字——对数边际可能性。

上面的推导，尤其是第(6)行，证实了我们不能计算这个 KL 散度。这是因为它要求我们计算对数边际可能性log p(y)，但是 p(y) 是一个难以处理的积分，它首先阻止我们使用贝叶斯规则来计算后验概率:

既然无法计算 logp(y) ，那就勇敢一点，把它去掉，尽量把剩下的两项最小化。

等价的，我们可以最大化的否定这两个项。我们给这种否定命名为证据下界，或 ELBO :

第(2)行否定了上面 KL 散度公式的前两项，去掉了第三项 log p(y) ，我们无法计算。

第(3)行和第(4)行使用 log 的属性和期望的线性来重新组织术语，即对于两个随机变量 a 和 b ，使用E【a+b】= E【a】+E【b】。

行(5)认识到第二项是变分分布 q(f) 和之间的 KL-散度，GP 先于 p(f) 。

第(6)行将第一项中的期望扩展到其实际的数学定义——积分。

从现在开始， ELBO 将是我们最大化的目标函数，以使变分分布 q(f) 尽可能接近后验 p(f|y) 。

为什么可以删除 log p(y)？

随着 EBLO 的定义，我们有了这个等式:

注意 ELBO 前面的减号。

我们的理想目标是最小化 KL(q(f)||p(f|y)) ，但是因为我们不能计算 log p(y) ，我们最小化 -ELBO ，或者等价地，最大化 ELBO 作为的替代*。*

你可能会问为什么可以去掉 log 边际可能性 log p(y) ？回答这个问题，我可以引用“所有的模型都是错的，但有些是有用的”。我们知道丢弃 log p(y) 是错误的，但是如果我们保留它，我们就无法计算公式。

*玩笑归玩笑，我认为你的问题的潜在含义是:丢弃 log p(y) 我们会失去什么信息？这种信息丢失是可以容忍的吗，这样剩下的模型仍然有能力发挥作用？为了回答这个新问题，让我们考虑一下 log p(y) 包括哪些信息——它是由先前的 p(f) 加权的可能性 *p(y|f】、之和:

下面再次显示的 ELBO 公式也包括出现在第一项中的可能性信息和出现在第二项中的先验信息:

因此，通过丢弃 log p(y) ，我们并没有失去一切，这是重要的一点。当然，剩余的信息没有以加权和的形式出现，所以我们有一些信息损失。这就是我们使用可计算的替代公式而不是理想的、但不可计算的公式所付出的代价。

ELBO 背后的直觉

从上面看第(6)行，我们对 ELBO 有一个直观的解释，复制到这里:

*第一项是相对于随机变量 f 的期望似然性，现在来自变分分布 *q(f)。第二项是 q(f) 与 GP 前 p(f) 之间的 KL 背离——注意，不是后 p(f|y) ，记住一个 KL 背离是非负 。

我们想要最大化 ELBO 。要做到这一点，我们需要期望的可能性大，KL 距离小。这个数学意义很好地符合我们的意图:

• 我们希望我们的模型 q(f) 能够通过很大的可能性很好地解释训练数据。
• 同时，我们希望 q(f) 通过 q(f) 和 p(f) 之间的小 KL 距离，保持接近我们的先验。这是因为我们相信我们先前的模型如何生成训练数据。后验仅重新加权来自前验的函数，使得很好地解释训练数据的函数获得更高的概率。理解 KL(q(f)||p(f))项的另一个角度是它是一个正则项。

ELBO 提到先前的p(f；l，σ ) 和变分分布q(f；μ, Σ).因此它仍然是这四个模型参数的函数。在第一项中，对数似然项只提到了变分参数 μ 和σ。第二项，即 KL(q||p) 项，提到了所有四个模型参数。**

你可能会想，既然我们放弃了我们无法计算的项——在KL(q(f)| | p(f | y)】公式中的 logp(y) 项，为什么我们不能在意识到有一个难以处理的积分的第一时间放弃它呢？那时我们试图用贝叶斯法则来计算后验概率:

并且我们意识到 p(y)=∫p(y|f)p(f)df ，也就是分母，很难处理。但是我们不能把它放在这里，因为我们会得到一个没有分母的分数，这是不良的。

我们现在可以去掉 logp(y)) 项，因为 KL-divergenceKL(q(f)| | p(f | y))由其他项组成。

你可能会奇怪，既然 KL 散度衡量的是两个分布之间的相似性，为什么我们选择做 KL(q(f)||p(f|y)) ，而不是 KL(p(f|y)||q(f)) ？

• 第一， KL 散度不对称，意思是 KL(q(f)||p(f)) 不等于 KL(p(f)||q(f)) 。
• 第二，如果我们使用 KL(p(f)||q(f)) ，我们会有问题:

从第(3)行我们可以看到，两个积分都提到了我们无法计算的后验 p(y|f) 。

好了，现在我们有了关于模型参数最大化的 ELBO 函数。因此 ELBO 是我们的优化目标。寻找具体值以最大化该目标的过程称为参数学习*。*

从对数边际似然角度看

您可能还记得，在高斯过程回归设置中，参数学习最大化的目标函数是 log 边际似然 logp(y) ，在该设置中， logp(y) 是易处理的。

现在在分类中，我们最大化 ELBO ，这两个目标函数是如何关联的？事实上，你可以从对数的边际可能性观点中推导出 ELBO :

线(1)是对数边际可能性的符号。

第(2)行将定义 p(y) 扩展为一个积分。

第(3)行引入了数量 q(f)/q(f) ，其计算结果为 1。所以这个量不会改变积分的结果。这是一个著名的操作——它将变分分布 q(f) 带入图像。

第(4)行把第(3)行中的积分改写成关于来自变分分布 q(f) 的随机变量 f 的一个期望。

第(5)行应用詹森不等式将 log 函数从预期之外移到预期之内。由此产生的期望是 ELBO 的公式。我们使用詹森不等式给我们计算带来方便——一个对数的期望比一个期望的对数更容易操作。

ELBO 是一个下界

请注意在应用詹森不等式后，从等号变为大号或等号。它揭示了 ELBO 小于我们真正希望最大化的 log 边际可能性 log(p(y)) 。这个边际可能性也被称为模型证据。而我们的 ELBO 比这个证据低，因此得名证据下界。

你可能会问，这个 ELBO 下界和我们真正想要最大化的 log 边际似然 log p(y) 有什么区别？我们已经从前面的等式中知道了这种差异:

差异是变分分布 q(f) 和真实后验分布 p(f|y) 之间的 KL 散度。重新排列一些术语后，我们有:

注意一个 KL 散度是非负量，所以 KL(q(f)||p(f|y)) ≥0 。这个等式证实了最大化 ELBO 等价于最小化 q(f) 和 p(f|y) 之间的 KL 散度，因为更大的 ELBO 意味着更小的 KL(q(f)||p(f|y)) 。为了正确理解这一点，我们需要了解一些微妙之处:

1. log(p(y)) 是提及核参数 lengthscale l 和信号方差 σ 的表达式。
2. ELBO 是提及核参数以及变分参数 μ 和*σ的表达式。*
3. KL(q(f)||p(f|y)) 是提及核参数以及变分参数 μ 和*σ的表达式。*

*通过改变 l、σ、μ、和σ的值来最大化 ELBO 将同时改变 log p(y) 的值。在参数学习过程中，不能将 log p(y) 视为常数。但是无论参数学习过程为 *l，σ，选择哪个值，一旦 ELBO 在这个模型参数值配置中最大化，这个配置就给你一个最小化的 KL(q(f)||p(f|y)) 。

ELBO 的解析表达式

我们再来看看爱尔博公式:

ELBO 是一个以模型参数为自变量的函数。我们需要推导出 ELBO 的解析表达式，因为我们想要使用梯度下降法来找到使 ELBO 最大化的那些模型参数的具体值。梯度下降仅适用于分析函数，因为它计算该函数相对于其参数的梯度。

KL 项已经是解析的了

由于 q(f) 和 p(f) 都是多元高斯概率密度函数，它们的 KL-散度的解析形式是已知的。这是选择 q(f) 做多元高斯的好处之一。稍后，我会给你看这个 KL 项的解析表达式。

可能性项不是分析性的

可能性项是一个积分，所以还不是分析性的。让我们稍微操作一下:

以上表达式由结构相同的 n 项组成。为了导出整个似然项的解析表达式，我们可以关注任意项，比如第 I 个项项。下面是一个非常重要的引出:

第(1)行是似然项中第个项的公式。它是对随机变量向量 f 的一个 n 维积分。 n 是训练数据点的数量。这是因为这个积分是计算【log(p(yᵢ|fᵢ】)相对于变分分布 q(f) 的期望，这是一个 n 维多元高斯。注)在 fᵢ 中只提到了单个随机变量，而 q(f) 提到了所有 n 个随机变量。

*第(2)行将随机变量向量 f 扩展成单个随机变量 f₁ 到 *fₙ.在 f 中有 n 个随机变量，所以有 n 个一维积分。

第(3)行引入一个名字 f₋ᵢ 来表示随机变量【f₁，f₂，…，fᵢ₋₁，fᵢ₊₁，…，fₙ】的向量。这是一个长度为 n-1 的向量。包含 f 中除 fᵢ.外的所有随机变量我们引入 f₋ᵢ 来简化公式。

第(4)行将 n 个积分重新分组为两个积分:一个在单变量 fᵢ 上，另一个在随机变量向量 f₋ᵢ 上。原因是我们想要应用多元高斯边缘化规则来计算 f₋ᵢ 上的内部积分。

第(5)行应用多元高斯边际化规则来计算对 f₋ᵢ 的内部积分，结果为 q(fᵢ).令人惊讶的是，我们知道 q(fᵢ的解析表达式。这一步需要更多的解释。

召回q(f；μ，σ)是均值向量 μ 和协方差矩阵σ的多元高斯分布。这里我显式地展示了 f 、 μ 和σ*的条目 :*

第(3)行的上述推导将分布 q 中的随机变量重新组织到 fᵢ、f₋ᵢ 块中。在这一行:

1. μᵢ 是一个浮点值；它是随机变量 fᵢ.的平均值
2. μ₋ᵢ 是一个长度为 n-1 的浮点向量；它是向量 f₋ᵢ.中随机变量的平均值
3. **σᵢ,ᵢ是一个浮点值；它是随机变量 fᵢ 的方差
4. **σᵢ,₋ᵢ和σ₋ᵢ,ᵢ是大小为 (n-1)×(n-1)的矩阵；它们是互相转置的。它们代表随机变量 f₋ᵢ 和随机变量 f₋ᵢ.之间的协方差
5. **σ₋ᵢ,₋ᵢ*是一个大小为 *(n-1)×(n-1)，的矩阵，它是 f₋ᵢ.中随机变量之间的协方差

随着分布 q 以这种方式重新表述，让我们回过头来看看随机变量 f₋ᵢ 的内部积分:

这个积分从联合分布 q(fᵢ，f₋ᵢ) 中整合出随机变量 f₋ᵢ ，得到 q(fᵢ) ，一个单一随机变量 fᵢ.的分布这就是概率论中边缘化的含义。

由于 q(fᵢ，f₋ᵢ) 是多元高斯分布，我们可以应用多元高斯边际化法则，繁琐地推导出 q(fᵢ) 的解析表达式，相当于计算上述积分。这是选择 q(f) 为多元高斯分布的另一个好处。

多元高斯边际化规则表示边际化分布 q(fᵢ) 仍然是高斯分布，其均值和方差等于联合分布中用红框突出显示的部分:

突出显示的部分是随机变量 fᵢ 的一维高斯分布，具有均值 μᵢ 和方差σᵢ,ᵢ。μᵢ* 是 μ 向量的第I项，σᵢ,ᵢ是**矩阵σ的对角线*的第I项。概率密度函数 q(fᵢ) 为:

由此，我们可以看出第 i 个似然项为:

这个函数在 μ 中只提到一个 μᵢ 条目，在σ中只提到一个σᵢ,ᵢ条目。推广到所有的似然项，我们知道总体似然提及向量 μ 中的所有条目，并且仅提及σ的对角线条目。**

似然项提到了哪些模型参数？

重要的是，在似然项中只提到了变分参数 μ 和变分参数σ中的对角线项。在似然项中没有提到σ中的非对角线条目。在似然项中也没有提到核参数 l 和 σ 。**

可能性是 n 个一维积分的总和

上述公式揭示了第个似然项是单个随机变量 fᵢ 的一维积分。推广到似然性中的任意项，我们知道它们都是一维积分，每个积分针对 f 中不同的随机变量。

因此，如果我们能够计算出代表第个*个可能性项的积分结果的解析表达式，我们就能够计算出整个可能性的解析表达式。KL 项已经是解析的，我们将得到整个 *ELBO 的解析表达式。这将使我们能够使用梯度下降进行参数学习。

第 i 个似然项是一维积分的事实非常重要。因为这个可能性项是一个积分。我们之前说过，它是一个运算，而不是一个我们可以求值的解析表达式。这种一维积分也很难计算。但是由于它是一维的，我们可以应用一种叫做高斯求积的方法，以解析的方式来近似这个积分的结果。

高斯正交

高斯求积法则，或者更具体地说，高斯-埃尔米特求积法则，如下所示，表示左手侧的特定结构的积分结果可以通过右手侧的 m 项的和来近似，其中 F(t) 是 t 的函数:

在右侧， wⱼ 和 tⱼ 是预先计算的，我们可以从高斯正交表中查找它们。 F(tⱼ) 是 F 在 tⱼ 位置的评价。右手边是一个解析表达式，因为我们可以查找 wⱼ和 tⱼ，我们可以计算 F(tⱼ).左侧和右侧都计算出一个数字。

请区分我们介绍的两种近似值:

1. 之前我们引入了变分分布 q(f) 来逼近真实的后验概率 p(f|y) 。我们迫使 q(f) 接近 p(f|y) 的方法就是最大化 ELBO 。这是一个分布到分布的近似值。
2. 现在我们需要推导出ELBO中第 I 个似然项的解析表达式。由于积分很难计算，我们想用高斯求积来近似这种积分的结果。这是一个数字对数字的近似值。

现在，请接受高斯正交近似是有意义的。事实上，这个近似值通常非常接近积分的真实结果。从这个意义上说，高斯正交近似是真实积分值的估计量。因为近似值与积分结果不同，所以这个估计量有一个偏差，由近似值与真实积分结果之间的差异来衡量。如果选择足够大的 m ，偏差可以任意小，代价是更多的计算。**

高斯求积背后的直觉

我会写另一篇文章来解释为什么高斯求积是有意义的。但是现在为了让你好奇的心平静下来，这里有一个简单的解释。

您可以将积分结果理解为无限多个等宽矩形的总和，如下图所示。每个矩形都有一个无限小的宽度，其高度是被积函数在该矩形 x 位置的值。

图片来自此处

我们可以通过在不同的 x 位置仔细选择几个(而不是无限多个)矩形，并给这些矩形不同的宽度，来近似积分的结果。在函数缓慢移动的区域，矩形较宽，在函数快速移动的区域，矩形较窄。下图说明了这一想法。

图片来自这里

高斯求积公式中的加权和结构正好实现了上述思想。求积理论告诉你如何选择位置 tⱼ 来计算被积函数，以及这个矩形的宽度 wⱼ 应该是多少。这些位置和宽度已经在的表格中计算出来了。你只需要决定你想要多少个矩形，那就是 m 的值。你想要的矩形越多，近似越精确，但是计算越昂贵。

很神奇，不是吗？

用高斯求积逼近第 I 个似然项

现在，让我们使用高斯求积来推导近似第个似然项的结果的解析表达式。

再看集成:

这不是高斯求积法则左侧尺寸所要求的格式:

我们需要将这个积分转换成高斯求积规则所要求的格式。很多时候，我们在数学中做的事情是模式匹配。

由于第 I 个*似然项中的积分是关于一个高斯随机变量 *fᵢ ~ 𝒩(μᵢ，σᵢ,ᵢ)*的一个期望，我们先应用积分中的变量变化规则，将其转化为关于一个标准高斯分布 *𝒩(0，1)的一个期望。我们引入一个新的一维随机变量 u 来自标准高斯分布: u~𝒩(0，1)。

我们可以将随机变量 fᵢ 定义为:

使用高斯线性变换规则，很容易验证以这种新方式定义的 fᵢ 仍然来自分布 𝒩(μᵢ，σᵢ,ᵢ):

你可以看我的另一篇文章揭开 Tensorflow 时间序列的神秘面纱:局部线性趋势了解更多关于高斯线性变换规则的信息。

利用这个新的随机变量 u ，我们将第个似然项转换为:

线(2)扩展了 fᵢ.的概率密度函数

第(3)行应用变量改变规则将积分变量从 fᵢ 改变为 u 。

第(4)和(5)行简化了公式。在第(5)行，出现了标准高斯分布的概率密度函数。

第(6)行认识到，现在第个似然项是关于标准高斯随机变量 u 的期望。

我们意识到第(5)行的公式仍然不是我们应用高斯求积所需的格式。但是理想的格式只是另一个变量的变化。我们引入一个新变量 t :

并将第 I 个似然项变换为:

我们最终将第一个似然项转换为高斯求积规则所要求的格式，其中:**

*在变量变化的这两个应用中，所有新旧变量 *fᵢ、*和 *t、的定义域都是从-∞到+∞，所以这些运算不会改变积分的边界，所以我们一起忽略了边界。

通过应用我们期待已久的高斯求积公式，我们可以得出:

第(1)和(2)行是高斯求积公式。

线路(3)插入 F(t) 。

第(4)行写出了所有的可能性。让我们仔细阅读它，因为它揭示了重要的信息:

1. 这个公式是解析的，因为我们可以在将 Yᵢ 插入 yᵢ 后评估其值，决定值 m ，并用高斯求积表中的值替换所有 wⱼ 和 tⱼ 。它计算出一个实数值。
2. 由于每个第 I 个似然项由 m 个项组成，所以整个似然项由 m×n 个项组成，其中 n 为训练数据点的个数。这个好像是很多项，比如说 m 是 40，训练点数是 1000，我们对于似然项就有 40000 项。但这是可以承受的，因为项数与 n 成线性比例关系，比例系数为 m。对于高斯过程来说，任何线性变化都是好消息。
3. 这个表达式是一个函数，它提到了变分参数 μ 和来自σ的对角元素。因为这个表达式是 ELBO 的一部分，ELBO 是我们要最大化的目标函数，它确实满足了我们对目标函数的要求——它至少需要提到我们的一些模型参数。 ELBO 中的似然项只提到了部分变分参数，并没有提到来自 GP 先验的核参数。请记住这一点，我们稍后会详细说明。

求积公式中的常数

为了让您对高斯求积表有所了解，当我们选择 m=10 时，位置 tⱼ 和权重 wⱼ 的值为:

使用的项数越多(m越大*)，近似值就越精确，但计算成本也就越高。经验法则是使用 m=40 。*

恭喜你！我们推导出了整个 ELBO 的解析表达式，我们为参数学习做好了准备。在此之前，一个负责任的机器学习实践者会问几个问题。

高斯求积是一个好的近似吗？

由于我们使用高斯求积为我们提供了一个解析表达式来近似第 I 个 T2 似然项中的积分结果，我们需要知道这种近似有多好，以及它的计算负担能力如何。

高斯求积非常精确。你可以看看视频高斯正交的魔力——比第二好的东西好十亿倍来说服自己。事实上，高斯求积的偏差可以忽略不计。高斯正交近似结果与积分的真实值之间几乎没有差别。

高斯求积是一个确定性的过程，在近似值中没有方差——您调用高斯求积方法任意次，都会得到相同的结果。

因此，高斯求积是第个似然项的一个很好的近似。

为什么所有的模型参数都可以从训练数据中学习？

我们看到第 I 个似然项的解析表达式只提到了变分参数 μ 和另一个变分参数σ的对角线上的条目。也没有提到内核参数 l 和 σ ，这两个参数来自 GP previous。**

此外，我们注意到在 ELBO 中，只有可能性项提到了训练数据 Y 。你可能会问，σ的非对角条目和内核参数 l 和 σ 如何从数据中学习？这是一个非常重要的问题。我们的训练数据与所有模型参数之间的联系是一个必要条件，使得这些模型参数可以通过梯度下降从训练数据中学习。

要回答这个问题，我们需要考察 ELBO 中另一项的解析表达式，变分分布q(f；μ，σ)和先验的p(f；l，σ )* 。*

首先要注意的是矩阵求逆 K⁻ ，其中 K 是来自 GP 先验的协方差矩阵。K 的形状为 n×n ，其中 n 为训练数据点的个数。所以这种倒置是非常昂贵的。这也意味着，我们的模型不能扩展到大型数据集。

线(2)给出了在q(f；μ，σ)和p(f；l，σ )* 。两者都是相同随机变量 f 的多元高斯分布。计算有点复杂，请接受这是可以做到的，而且是分析性的。 n 是训练数据点的数量。 tr 代表跟踪算子，它是一个方阵对角线上的元素之和。 K 是 GP prior 中定义的 n×n 协方差矩阵。 det(K) 代表 K 的行列式。 K⁻ 是 K 的逆。而 m(X) 是我们给 GP 先验的均值函数，通常情况下， m(X)=0 。*

第(3)行将零均值函数 0 插入公式以简化它。

让我们用两个训练数据点 (X₁，Y₁)，(X₂，Y₂) 来研究第(3)行的公式。例如:

协方差矩阵 K 在训练数据中只提到了 X ，没有提到 Y 。作为协方差矩阵的 K 是对称的。并且 K 提到了内核参数 l 和 σ。

让我们看看我们的变分分布(这里，为了节省空间，我把协方差矩阵的一个条目写成σᵢⱼ*，而不是我以前使用的符号σᵢ,ⱼ😗

*有了上面的定义，我们可以使用微积分中的规则 : 计算行列式 det(K) ，*det(σ)和逆 K⁻

**

现在，我们可以写下 KL 散度公式中的各项:

我们只需要浏览一下这些项，就可以看到在 KL 散度项中:

1. 所有模型参数，包括来自变分分布的 μ 中的所有条目和σ中的所有条目(对角线以及非对角线)，以及来自 GP prior 的 l 和 σ ，都纠缠在一起。
2. 由于这种纠缠，似然项中提到的模型参数部分(即 μ 和σ的对角线条目)会影响似然项中没有提到的其余模型参数(GP prior 中的σ和 l 和 σ 的非对角线条目)。
3. 并且，训练数据 Y 可以影响所有模型参数。

这就是为什么所有模型参数都可以从训练数据中学习到。

模型可扩展性

从 q(f) 到 p(f) 的 KL 散度中的矩阵求逆 K⁻ 计算起来很昂贵。运算的次数，比如矩阵中的倍数或者两个数相加，与 n 成比例。 n 是训练数据点的数量。对于现代标准的微小数据集 n=1000 ，反演需要数十亿次运算。

*这意味着我们的变分高斯过程无法扩展到大型数据集——计算 *K⁻需要太长时间。我的下一篇文章将讨论稀疏和变分高斯过程，这是一个可以扩展到大型数据集的模型。

我们可以用高斯求积来计算后验概率吗？

你可能会问，我们在 ELBO 中用高斯求积来近似似然项，是否可以用它来近似贝叶斯法则中的边际似然？

变分分布 q(f) 和 GP 先验都是多元高斯分布。所以结构上，这两个整合是一样的。为什么不用高斯求积来近似边际似然呢？我完全理解这个问题——它也让我抓狂。

但是在 ELBO 似然项中的log(p(y | f)】和边际似然项中的 p(f|y) 有着至关重要的区别。边际可能性前面没有日志*。所以我们不能使用属性 log(ab) = log(a) + log(b) 将边际似然分解成更小项的和。没有这一点，我们不能将多元高斯边际化规则应用于边际似然，以将对 f 的 n 维积分减少到对单个随机变量 fᵢ.的单个积分之和*

因此，边际似然仍然是一个 n- 维积分，其高斯正交近似需要 mⁿ 项，即mn 项的幂——总项数随着训练数据点的数量异常增长。即使对于一个小的数据集来说，这也是一个不切实际的大数字。这就是我们不能在直接计算后验概率时使用高斯求积的原因。

本文的附录向您展示了为什么高斯求积需要 mⁿ 项来逼近边际似然。

参数学习

*现在我们已经推导出了 ELBO T1 的解析表达式。这个表达式是一个函数，其中我们的所有模型参数都是未知的——核心参数 *l，σ，和变分参数 𝜇，σ。参数学习使用梯度下降来最大化 ELBO 以找到这些模型参数的最佳值。这是梯度下降的直接应用。但是仔细思考这种优化提出了一些有趣的问题。

先总结一下我们的模式。

完整的模型

我们的二元分类模型包含以下组件:

1. 多元高斯先验p(f；l，σ )* ，其中引入了潜在随机变量向量 f ，以及核参数 l 和 σ。*
2. 伯努利似然 p(y|f) ，引入观测随机变量向量 y ，不引入任何参数。
3. 变分分布q(f；μ，σ)，其中引入了变分参数 μ 和σ。我们要求q(f；μ，σ)向后方靠近p(y | f；l，σ )* 在 KL-散度意义上。*

在数学语言中，模型显示为(红色的模型参数):

第(1)行是 GP 优先。第(2)行是可能性。第(3)行定义了由 Bayes 规则计算的真实但难以处理的后验概率。第(4)行引入了近似真实后验概率的变分分布。第(5)行使用 KL 散度来定义通过变分分布来近似真实后验概率的含义。这种 KL 差异也是难以处理的，所以我们最终最大化 ELBO 作为最小化这种 KL 的代理。

我们需要学习多少标量？

我们的模型有四个不同形状的参数(形状在张量形状的意义上):

• 长度刻度 l 为标量。
• 信号方差 σ 也是标量。**
• **q(f)的均值向量 μ 是一个长度为 n 的向量，所以它有 n 个标量项。
• **q(f)的协方差矩阵σ*是一个大小为 *n×n，的矩阵，它是对称的，所以有 n(n+1) 自由标量项。

这个模型总共有 2+n+ n(n+1) 个标量，这些标量的值需要通过参数学习来决定。不要忽略列出所有模型参数的任务——这是一个健全性检查，以便您理解模型。

我知道你想喊出来: 2+n+ n(n+1) 是很多标量，它比训练数据点的数量 n ！

既然大部分标量来自于变分分布 q 的参数，那我们能不能用更少的标量参数化 q ？是的，这就是平均场参数化。

q 的平均场参数化

我们的变分分布q(f；μ，σ)是完全多元高斯分布，因为σ是没有零点的完全协方差矩阵。如果我们要求σ是对角矩阵，那么 q(f) 就变成了平均场分布。

使用平均场 q(f) 的动机是有较少的变化参数要学习。代替来自全协方差矩阵的 n(n+1) 自由标量，我们现在只有来自平均场协方差矩阵对角线的 n 自由标量，因为它是对角方阵。在平均场配置中，我们需要学习的标量总数是 2+2n 。模型参数越少，参数学习过程越容易，这意味着优化停止得越快。

这种平均场近似的缺点是，它只能模拟更小的分布子集——约束所有随机变量相互独立的子集。然而，在现实中，真实后验概率中的随机变量通常是相互关联的。这意味着平均场近似与具有全协方差矩阵的变分分布相比是更粗略的近似。

这是我们必须决定的模型容量和计算成本之间的权衡。通常的做法是从平均场近似开始。如果学习模型在预测中表现不佳，请使用更昂贵的全协方差近似值。

当变分分布q(f；μ，σ)具有平均场结构，即σ是一个对角矩阵，在 f 中的所有单个随机变量都是相互独立的——它们之间没有任何关系。因此，关于一个个体随机变量的信息，比如说 fᵢ的和 fⱼ的不能用来推断另一个个体随机变量。你可能想知道我们的模型是如何做出预测的。

我们的模型可以在测试位置 X_ 进行预测，因为即使位于训练位置 X、的 f 中的随机变量之间的相关性为零，但是 f 与测试位置 X_ 的随机变量 f_ 之间的相关性不为零。这个非零相关性是我们通过先验强加的条件概率 p(f_|f)* 。它使模型能够通过使用来自 f 的信息来推断 f_* 的值。我们将在做出预测部分更清楚地看到这一点。*

参数太多的模型？

在统计了需要学习的标量数量后，我们意识到一个问题——我们的模型要学习的标量( 2+n+ n(n+1) 在全高斯配置中， 2+2n 在平均场配置中)比训练数据点的数量多，也就是 n 。

这很糟糕，因为模型现在有能力精确地拟合训练数据-模型可以只使用其参数来记住训练数据在哪里。但是这个模型不能推广到新的数据。人们用术语过度拟合来指这种模型任意地很好地拟合了训练数据，但是对测试数据做出了糟糕的预测的情况。

顺便说一下，模型具有足够的参数来任意地很好地拟合训练数据，并不意味着梯度下降可以找到那些参数的具体值来很好地拟合训练数据。模型的参数越多，梯度下降在参数值空间中导航就越困难。梯度下降可能收敛到一些甚至不适合训练数据的模型参数值。

理想情况下，我们希望模型有足够数量的参数和足够的模型结构(核函数、似然函数等)。)使得模型可以很好地总结训练数据的潜在特征，而不是记住训练数据在哪里。这种总结还使模型能够归纳出新的测试数据。但是从设计上来说，我们的二元分类模型，参数太多，不是一个好的模型。真扫兴。

事实上，你学习这个模型的真正原因是为了让你熟悉如何将变分推断应用于高斯过程模型。这是为了让你为实际的实用模型——稀疏变分高斯过程模型(SVGP 模型)做好准备。有关 SVGP 模型的更多信息，请访问:稀疏和变分高斯过程——当数据很大时该怎么办。现在我建议你先完成这篇文章。

没有随机梯度下降？

ELBO 由似然项和 KL 项组成。这两个术语都需要完整的数据来评估。如下所示，可能性项需要完整的 Y 。

而 KL 项，如下所示，需要所有的训练数据来计算:

第(3)行的公式提到了先前的协方差矩阵 K=k(X，X) 。 K 需要整个 X 来计算。

因此，在每个优化步骤中，梯度下降算法需要检查所有训练数据，以计算 ELBO 及其梯度。这就是为什么我们不能在我们的变分模型上使用随机梯度下降。随机梯度下降在每个优化步骤中只使用一批训练数据，因此它需要的内存少得多，而且通常收敛得更快。

好消息是，我将在下一篇文章中讨论的稀疏和变分高斯过程模型可以利用随机梯度下降。所以请继续关注，一切终将解决。

梯度下降压力太大？

我们使用梯度下降来最大化目标函数 ELBO 。如果我们仔细想想， ELBO 实际上在一个公式中编码了两种优化:

1. 找到核参数 l 和 σ 的值，使得真实后验很好地解释了训练数据。
2. 找到变分参数 μ，σ的值，使得变分分布q(f；μ，σ)非常接近真实的后验。

这对于梯度下降来说要求太高了。这在数学上意味着: ELBO 可以是高度非凹的，有很多局部极大值。梯度下降法作为一种局部优化方法，会卡在一个局部最大值上，不一定是一个好的局部最大值。那么我们如何优化这样的模型就成了它自己的艺术/科学:

• 使用哪些优化器？
• 我们如何降低学习率？
• 我们如何在优化开始前初始化模型参数？
• 我们如何处理数值不稳定性？对于高斯过程模型，当我们对一个矩阵求逆时，这意味着乔莱斯基分解失败。

这是一个很大的话题，我将在不久的将来发表几篇文章来解释它。

做预测

现在我们有了随机变量 f 的后验 p(f|y ) ，它在训练位置 X 。其实我们没有 p(f|y)。我们有类似于后验概率的变分分布 q(f) 。不过还是继续用 p(f|y) 吧，好像我们知道一样。

使用后验概率进行预测

先说清楚为什么需要后验来做预测。

我们的完整模型由三个随机变量组成， f_，f ，和 y 。 f_* 和 f 是潜在随机变量，它们没有相应的观测值。 y 是一个观察型随机变量，它有来自我们训练数据的观察值 Y 。*

通过 GP prior p(f_，f) 我们对随机变量 f 和 f_* 施加了一些特征。通过似然 p(y|f) ，我们建立了潜在随机变量 f 和观察随机变量 y 之间的联系。*

贝叶斯规则使用观测值 Y 以后验分布 p(f_，f|y) 的形式给出 f 和 f_* 的更新特征。这就是我们要使用后验概率进行预测的原因，后验概率会对先验概率中的所有函数进行重新加权，以便更好地解释训练数据的函数获得更高的概率。*

后验 p(f_，f|y) 是随机变量 f_* 和 f 的分布。但是我们只需要对 f_* 进行预测，换句话说，我们只需要对 f_* 进行分布。因此，我们对来自 p(f_，f|y)* 的随机变量 f 进行积分，得到预测分布 p(f_|y)* 。*

注意，即使我们的模型由三个随机变量 f_ 、 f 和 y 组成。我们不需要知道它们三个的联合概率密度，记为 p(f_，f，y)* 。在整篇文章中，这个联合分布 p(f_，f，y)* 从未出现过*。我们可以计算出我们需要的，也就是预测分布 p(f_|y) ，而不需要 p(f_，f，y)。这是一件好事，因为 p(f_，f，y) = p(f_，f|y) p(y)* 我们知道边际可能性 p(y) 不是一个容易处理的量。我们以定义的方式定义了我们的模型，以确保我们不必触及联合分布 p(f_，f，y)* 。*

高斯分布形式的后验

为了在测试位置 X_ 进行预测，我们需要推导随机变量 f_* 在测试位置的后验 p(f_|y)* 。*

第(2)行逆向使用概率论中的链式法则来拆分联合分布。

线(3)在 p(f_|f，y)的条件下下降 y 。我们可以这样做，因为 GP prior 只提到了 f 和 f_ 。所以给定 f ， f_ 独立于 y 。 y 在可能性中被介绍。*

线(4)用我们知道的变分分布 q(f) 代替我们不知道的后验 p(f|y) 。这是重要的一步:

• *q(f) 与后验 *p(f|y)具有相同的 API。还有，
• 在最小 KL 散度意义上，它的行为类似于后一种情况。

所以我们可以用 q(f) 代替 p(f|y) ，这就是变分推理的全部意义。**

在第(5)行，在积分符号内， p(f_)|f) 是通过对 GP 先验应用多变量高斯条件规则(此处为更多该规则的细节)而导出的条件分布:*

并且我们引入名称 A 和 B 来缩短公式，用:

A 和 B 都提到了训练位置 X 和测试位置 X_ 。它们建立了从 f 到 f_* 的连接。*

第(6)行应用多元高斯线性变换规则，通过仅使用来自q(f；μ，σ)不提 f 。

第(7)行将 f_ 的概率密度函数移到积分之外，因为它没有提到积分变量 f ，所以它可以被视为积分的常数。*

第(8)行将积分计算为 1，因为随机变量定义域上的概率密度函数的积分计算为 1。

结果是随机变量 f_ 在测试位置 X_的后验分布 p(f_|y) 。我们也称之为预测分布。如果 f_ 是一维的就是一维高斯分布，如果 f_ 是多维的就是多元高斯分布。*

预测分布

扩展了 A 和 B 的预测分布的完整公式为:

使用:

需要注意几件重要的事情:

首先，预测分布没有提到训练数据 Y 。您可能还记得在高斯过程回归模型(计算后验部分)中，预测分布的平均值是训练数据 Y. 的加权和，但是在我们的变分模型中，预测分布的平均值是来自变分分布 q 的平均值 μ 的加权和。在参数学习过程中， Y 的信息被完全吸收到 μ 的值中。

第二，预测分布提到了训练数据 X — 你可以看到 X 遍布在预测分布的均值和方差中。这是因为模型需要使用训练位置 X 和测试位置 X_之间的距离来计算变量 f(X) 和 f(X_) — 之间的相关性 X 和 X_ 之间的距离越远，相关性越小，模型给予 f(X) 的权重越小，以预测 f(X_)* 的值*

最后，我们理解我们的模型如何能够做出预测，即使变分分布 q 具有平均场结构，在这种情况下，σ是对角矩阵。我们的模型使用相关性矩阵 k(X_，X) 将训练位置 X 处的随机变量 f 和测试位置 X_* 处的随机变量 f_* 联系在一起。通过条件概率 p(f_|f)* 引入相关性 k(X_，X)；它不关心 f 中的单个随机变量是否相互独立。

伯努利分布形式的后验

在二元分类设置中，我们需要以伯努利分布的形式进行预测，而不是以高斯分布的形式。

为了构造伯努利后验概率，我们需要决定它的单个参数的值:

假设只有一个测试点，那么 X_ 是一个标量，而不是一个矢量。所以 f_* 是来自我们后验 p(f_|y)的一维高斯随机变量。以这种方式定义， g(f_) 是一个随机变量，但是伯努利分布需要一个标量作为其参数，因此我们计算 g(f_)* 相对于 f_😗 的期望值*

你会很高兴地意识到线(3)是一个积分，你可以用高斯求积来近似。

计算出 g_ 后， y_* 的伯努利概率密度为:*

或者相当于:

上述伯努利概率密度函数是我们对测试位置的最终预测。你可能会问，难道不是所有来自高斯过程模型的预测都带有不确定性度量吗？是的，伯努利分布确实有它的不确定性度量，继承自高斯后验分布。由于伯努利分布是离散分布，其不确定性被编码在概率参数 g_: 中*

• g_ 接近 1:模型确定事件 y_=1。**
• g_ 接近 0:模型确定事件 y_=0* 。*
• g_ 接近 0.5:模型不确定要预测哪个标签。*

记住 g_ 是一个导出量。它从后验的 p(f_|y)中得到它的值。**p(f _ * | y)*的均值和方差是决定 *g_值的两个东西(连同测试位置 X_ )。因此，它们决定了伯努利预测的不确定性。

实验结果

我用这个代码创建了一个变分高斯过程模型。它使用了 gpflow ，这是一个用于高斯过程模型的 Python 库。

下图显示了模型在参数学习后对我们的训练数据 (X，Y) 的预测。

*顶部子图显示了来自高斯后验 f(x) = *p(f|y)的预测。预测的形式是均值和方差。所以上面的子图用黑色显示平均值，用青色显示 95%的置信区间。可以看到，对于标签为 1 的数据点， f(x) 值在 4 左右，对于标签为 0 的数据点， f(x) 值在-4 左右。

中间的子图用蓝色显示挤压函数 g(x) 。我们可以看到，对于标签为 1 的数据点，它的值接近 1，对于标签为 0 的数据点，它的值接近 0。

底部的子图用红色显示了我们的训练数据。

这个情节和本文开头的手绘情节很像。

结论

这里我总结了我们的变分高斯过程模型的主要步骤:

1. 该模型使用多元高斯先验 p(f) ，并使用伯努利似然 p(y|f) 对二元观察值进行建模。
2. 当可能性是非高斯的时候，贝叶斯规则会导致一个难以处理的后验 p(f|y) 。
3. 使用变分分布 q(f) 来近似后验 p(f|y) 。我们通过要求 KL 散度 KL(q(f)||p(f|y)) 最小来定义一个好的近似。
4. *KL(q(f)||p(f|y)) 的公式也很难处理。所以我们使用 ELBO 公式作为 KL(q(f)||p(f|y)) 最小化任务*的代理。我们需要最大化 ELBO。艾尔博的公式很容易理解。我们需要推导出 ELBO 的解析表达式，这是一个提到所有模型参数的表达式。这样，我们可以使用梯度下降来最大化关于模型参数的 ELBO 以找到那些模型参数的具体值。
5. ELBO 由两项组成，似然项和 KL 项 KL(q(f)||p(f)) 。KL 项已经是解析的了。但是可能性项是一个很难计算的积分。
6. 似然项是一个 n 维积分，可以分解成 n 维一维积分。高斯求积为每个一维积分的结果提供了一个解析近似。
7. 梯度下降使 ELBO 的解析表达式最大化，以找到具体的模型参数值。
8. 最后，我们指出，我们的模型包含的参数比训练数据点的数量多。这不是一个实用的模型。它为我们接下来将要学习的真正实用模型——稀疏和变分高斯过程模型做准备。

支持我

如果你喜欢我的故事，如果你考虑支持我，通过这个链接成为灵媒会员，我将不胜感激:https://jasonweiyi.medium.com/membership。

我会继续写这些故事。

参考

变分高斯过程

附录:可以用高斯求积来近似边际似然吗？

在这里，我给出的推导表明，贝叶斯规则的边际似然的高斯正交逼近需要 mⁿ 项。目标是让你相信使用高斯求积来近似边际似然项是不实际的。我建议你仔细阅读这些推导，因为它们是非常好的贝叶斯统计实践。

这又是贝叶斯法则:

我们之前说过，我们不能用这个规则来计算后验概率，因为分母的积分是很难处理的。现在让我们尝试使用高斯求积来近似这个积分。

可能性 p(y|f) 是 n 伯努利概率密度函数的乘积，如我们之前推导并在此再次显示的:

由于在 p(y|f) 之前没有 log ，我们不能使用 log 的属性将似然分解为n**logBernoulli 概率密度函数之和，每个函数针对单个数据点。这将导致高斯求积规则产生指数数量的项。要理解为什么，让我们继续分析。

让我们通过应用变量变化规则来操纵这种集成:

第(1)行是我们想要导出解析表达式的分母中的积分。被积函数是伯努利似然 p(y|f) 乘以 GP 先验 p(f) 。

第(2)行写出了 GP 先验的分布，它是一个均值为 0、协方差矩阵为 K 的多元高斯分布。

第(3)行通过引入新变量 u 并将原始积分变量 f 重写为 Lu ，其中 L 是协方差矩阵 K 的乔莱斯基分解的下三角矩阵，从而在积分中应用变量变化规则。并且 u 是具有 0 均值和恒等协方差矩阵 1 的 n 维标准高斯分布。 0 是长度为 n 的 0 向量， 1 是形状为 n×n 的单位矩阵。

*第(4)行简化了第(3)行。我们需要更多关于如何导出这种简化的解释，但首先，我们验证重新定义为 Lu 的 f 仍然来自同一个高斯分布 *𝒩(0，k)。由于 f 是从 u 的线性变换，我们用高斯线性变换规则把 f 的分布写成从 u 的分布，这是一个多维的标准正态分布 𝒩( 0 、 1 。

有关高斯线性变换规则的更多详情，请参见揭秘张量流时间序列:局部线性趋势，搜索“线性变换”。

现在我们已经验证了 f 仍然来自 GP prior。为了说明我们如何在上面的第(4)行得到积分，让我们用一个有 3 个数据点的例子。因此

扩展到:

为了应用变量变化积分规则，我们需要计算雅可比矩阵，然后计算雅可比矩阵的行列式:

*以上所有推导都很简单。在第(4)行引入名称 a 、 *b、和 c 以缩短公式，因此所有内容都在一行中。

这个推导表明，雅可比矩阵等于来自乔莱斯基分解的下三角矩阵 L 。所以 det(J(u)) = det(L) 。

现在让我们看看变量变化 f=Lu 对 GP 优先 p(f) 的影响:

第(2)行使用属性 det(K) = det(L)。这是来自乔莱斯基分解的一处不错的房产。

线路(3)插入 f = Lu 并使用同样来自乔莱斯基分解的属性:

第(4)行简化了公式。它将 det(L) 移出原来的分母，加上det(1)*，等于 1 *。这是因为我们要构造单位多元高斯的公式。

第(5)行认为公式等于概率密度 u ，由 1/det(L) 缩放。而我们将愉快地认识到 1/det(L) 取消了来自det(J(u)】的 det(L) 。

在我们结束变量变化推导之前，我们需要验证变量变化运算前后的积分边界。

原始积分将 f 作为其积分变量。 f 是一个高斯随机变量的向量。所以 f 中的每个随机变量都有一个全实直线的定义域。我们应用变量变化规律后，积分变量变成 u ，也是一个多变量高斯随机变量，全实域。这就是为什么我们对导数中的边界是隐式的。

使用高斯求积

现在我们有了以下积分，我们希望使用高斯求积来给出近似其结果的解析表达式:

第(1)行是要近似的积分。第(2)行引入了 name h(u) = p(y|Lu) 来缩短下面推导的公式。

在我们的 3 个数据点示例中， u 是一个三维标准正态分布，因此上面是一个三维积分。但是高斯求积法则只适用于一维积分。

标准多元高斯 𝒩(u 的好处是: 0 ， 1 ) 是所有的个体随机变量 u₁ ， u₂ ， u₃ 因为恒等式协方差矩阵 1 而相互独立。于是联合分发𝒩(u；* 0 ， 1 ) 分解成:*

通过这种分解，我们有:

第(2)行将 u 扩展为单个随机变量及其概率密度函数。可能性 h( u₁ 、 u₂ 、 u₃) 是所有三个随机变量的函数。

第(3)行将关于随机变量 u₁ 的积分分成一对括号。因为我们将首先解决这个最内部的整合。

这就是:

第(2)行扩展了𝒩(u₁的概率密度函数；0，1) 。

第(3)行重新组织指数函数中的项，为下一行中变量的变化做准备。

第(4)行应用积分中的变量变化规则，新变量 t₁ 使得 u₁=√2 t₁.

第(5)行将公式简化为我们可以应用高斯求积规则的形式。

第(6)行应用高斯求积规则来近似具有 m 项的积分。

术语数量呈指数增长的原因

注意第(6)行，即使随机变量 u₁ 已经被高斯正交“积分”出来，函数 h 仍然包含随机变量 u₂ 和 u₃.这就是指数数量的术语的来源。因为现在我们必须应用高斯求积法则，一次又一次，一个接一个地积分随机变量。每次，我们积分出一个随机变量，接下来是 u₂的和 u₃的。每一次，被积函数变得越来越长。**

例如，对于相对于 u₂ 的集成，我们有:

由于双重求和，上述表达式有 m 项。

我们需要应用高斯求积规则三次来整合所有三个随机变量。之后，我们得到了整个边际似然项的最终解析表达式:

这个最终近似值有 m 项。

为了将其推广到具有 n 个训练点的数据集，我们将在近似边际似然项的分析表达式中有个 mⁿ* 个求积项。因此，术语的数量与数据点的数量成指数关系。它很快变得大得不切实际。这就是为什么我们不能使用高斯求积来近似贝叶斯规则中的边际似然项。*

贝叶斯多元高斯混合模型中的变分推断

原文：https://towardsdatascience.com/variational-inference-in-bayesian-multivariate-gaussian-mixture-model-41c8cc4d82d7?source=collection_archive---------16-----------------------

变分推理是贝叶斯统计中的一种近似推理方法。给定一个模型，我们经常想要推断它的后验密度。然而，关于后验密度的精确推断需要通过所有不相关的参数进行边缘化，这通常在计算上是难以处理的。VI 的主要思想是使用简单分布来近似查询密度。为了简化问题，我们用一个分布来近似每个单独的目标参数，并假设我们用于近似的分布是独立的。这种方法叫做平均场近似法。

等式 1

说到近似，我们自然会想到 KL 散度，它衡量两个分布之间的相似性。但是 KL 散度本身就包含了我们要推断的后验密度。因此，这又是一个棘手的问题。然后我们介绍**证据下界(ELBO)。**可以证明 KL 散度的最小化等价于 ELBO 的最大化。

等式 2

现在优化是可能的。我们所做的是，我们遵循与期望最大化相似的程序( EM )，一次更新一个近似分布。原来更新函数的表达式相当简单，如下图所示。

等式 3

注意，z 是我们想要近似其分布的参数。

上面提到的方程是我们进行变分推理所需要的。让我们看看它们在贝叶斯多元高斯混合模型中的具体应用。该模型的图示如下所示。

每个节点的分布如下:

这个模型的直觉是，观察值或数据是从高斯分布生成的。这个高斯分布的参数由另外两个分布控制。首先，我们有 k 个μ选择，它也是由高斯分布产生的。其次，我们有一个分类分布，它决定 k 个生成的μ中的哪个均值负责生成数据。概括地说，生成数据的分布是多元的，这意味着我们的观察值 x 不止一个维度。我们用 p 来表示数据的维数。

这里的异议是，我们想用平均场近似来近似 p(μ，c)。

注意，根据概率图形理论，这里的μ和 c 不是独立的，因为数据是观察到的。然而，当我们设计近似函数时，我们可以使我们的近似分布相互独立。如你所见，这里的第一件事是找到所有 q 分布的形式。虽然我们可以使用任何我们想要的发行版，但是我们希望尽可能地简化它们。自然，我们希望使用与原始参数相同的分布(实际上，为了数学上的方便，我们希望先验与后验共轭。但这里不讨论。).所以，我们把近似分布设计成，

因此，我们的工作变成了优化 mk、sk 和ϕi.

让我们看看我们的 ELBO 长什么样

对于 ELBO 中的每一项，我们可以导出一个封闭形式。

第一期

因为我们假设，

所以，

第二学期

第三期

因为

所以，

第四项

期限 5

对 ELBO 的上述计算表明它是易处理的，并且可以通过监视它是否在增加来跟踪参数的更新。

现在我们使用等式 3 来导出参数的更新函数。

对于 q© ，

通过使用来自图形模型的联合分布，

所以，

对于 q(μ) ，

而且，

所以，

因为

所以，

还有，

然后，我们有了所有近似参数的更新函数。

我们可以在数据维数为 2 的情况下测试我们的理论。结果如下所示。

vcpkg:用于机器学习的 C++包管理器

原文：https://towardsdatascience.com/vcpkg-c-package-manager-for-machine-learning-16c3b18c1ef2?source=collection_archive---------27-----------------------

使用 vcpkg 安装和使用 C++ ML 库

使用 TensorFlow C API 的控制台应用程序

Vcpkg 是微软提供的一个很棒的 C++包管理器。它极大地简化了 Windows、Linux 和 macOS 上第三方库的获取和安装。

要在您的 C++项目中使用它，您首先需要将存储库克隆到您的本地机器上。仓库位于这里。克隆时，选择一个在路径上容易找到的文件夹很重要。类似于 C:\dev\vcpkg 的东西。在命令提示符下，您可以这样做:

git clone [https://github.com/microsoft/vcpkg.git](https://github.com/microsoft/vcpkg.git)

将存储库克隆到本地文件夹后，转到该文件夹并发出以下命令来构建可执行文件(vcpkg.exe):

.\bootstrap-vcpkg.bat

现在您已经构建了可执行文件，您可以使用它将 vcpkg 集成到您的构建环境中(如 Visual Studio 2019)。从命令行执行以下操作:

.\vcpkg integrate install

假设您想在您的计算机上安装 OpenCV4，这样您就可以将它集成到您的 C++程序中。你会这样做:

.\vcpkg install opencv4 --triplet x64-windows

三元组描述了您希望为其安装软件包的架构和环境。如果这是您安装的第一个软件包，vcpkg 将下载并安装一个兼容的 cmake，然后使用它来构建软件包。一旦完成，它将计算出 opencv4 的所有依赖项，并首先构建它们。最后，它将构建 opencv4 包并将其存储在包缓存中。一旦你安装了一个包，你就可以在你的 c++项目中使用它。

创建自己的包

创建自己的包非常简单。我在这里举一个例子，我用 TensorFlow 团队创建的预打包文件为 TensforFlow C API(见 https://www.tensorflow.org/install/lang_c)创建了一个。在命令提示符下，执行以下操作:

.\vcpkg create tensorflow-c-gpu https://storage.googleapis.com/tensorflow/libtensorflow/libtensorflow-gpu-windows-x86_64-2.3.0.zip libtensorflow-gpu-windows-x86_64-2.3.0.zip

在这里，您已经将您的包命名为 tensorflow-c-gpu ，并告诉它从哪里下载 zip 文件，以及您将在本地将其命名为什么。发出这个命令将创建一个 cmake 文件，用于构建您的包。在我们的示例中，该文件将是ports \ tensor flow-c-GPU \ tensor flow-c-GPU . cmake。每个包还需要一个名为 vcpkg.json 的元数据文件。要开始编辑这些文件，您可以从命令提示符发出以下命令:

.\vcpkg edit tensorflow-c-gpu

对于我们的例子来说, vcpkg.json 文件应该是这样的:

您还应该编辑创建的 cmake 文件，以使用下载的 zip 文件，并将它放在正确的位置。下面是我的 tensorflow-c-gpu.cmake 文件的样子:

从这个文件中可以看出，这个包只支持 Windows 上的 x64 版本。但这就是我现在所需要的🙂

使用包

让我们创建一个使用这个包的控制台应用程序。在 Visual Studio 中创建一个 C++控制台应用程序。我把我的叫做 TfConsole 。我将在本文末尾包含一个源代码链接。

一旦你创建了这个应用程序，添加一个名为 vcpkg.json 的新文件。这将描述应用程序及其所有依赖项:

你可以看到我们已经指定这个项目需要我们之前定义和安装的 tensorflow-c-gpu 项目。

您需要更改项目的构建属性，以包含此清单文件:

TfConsole 生成属性

将使用 Vcpkg 清单选项更改为是，如此处所示。

如果您在配置属性中没有看到 Vcpkg 选项，这意味着您还没有运行前面显示的“integrate install”命令。

您还需要添加 tensorflow.lib 库作为链接器的输入:

链接器属性

在 main.cpp 中，我们将对 tensorflow C API 进行简单的包含，然后使用它来显示 tensorflow 版本号:

需要注意的一点是，该包包含 TensorFlow 2.3.0 的 GPU 版本，因此您必须在您的机器上安装 CUDA 10.1。

如果一切顺利，您将在控制台输出中看到版本号。应该是 2.3.0。

控制台应用程序正在运行

结论

我认为微软的 vcpkg 管理器是一个很棒的包管理器，越来越多的包被添加到它里面。我还向您展示了如何添加您自己的，一旦您知道自己在做什么，这可能会非常简单。我正计划写一系列关于在 Windows 环境下生产中使用机器学习模型的文章。它将涵盖一个 Windows 前端(这将是一个 WPF C#桌面应用程序)使用各种不同的后端。其中一个后端将使用 tensor flow C API——由于没有相关的文档和很少的使用示例，它被证明是非常核心的！但是它也非常快(它的推理速度比我的 Python 代码还要快)。

观看此空间…

资源

1. https://docs.microsoft.com/en-us/cpp/build/vcpkg?view=vs-2019
2. VCP kg GitHub:https://github.com/Microsoft/vcpkg
3. TF console:https://github.com/ianormy/TfConsole

预测时间序列的向量自回归

原文：https://towardsdatascience.com/vector-autoregressive-for-forecasting-time-series-a60e6f168c70?source=collection_archive---------10-----------------------

VAR 捕捉发展和相互依赖性

包含多个时间序列的经济计量模型

作者图片

https://sarit-maitra.medium.com/membership

V 向量自回归(VAR)时间序列模型在计量经济预测模型中有着广泛的应用；VAR 可以捕捉多个时间序列之间的演化和相互依赖关系。VAR 中的所有变量都是对称处理的，每个变量都包含一个方程，根据其自身的滞后和模型中所有其他变量的滞后来解释其演变。我们可以称之为交易策略的科学方法。

使用科学方法进行交易策略设计的优势在于，如果策略在前一段盈利期后失败，则有可能重新审视最初的假设并重新评估它，从而潜在地产生新的假设，使策略重新盈利。

布伦特原油有四个时间序列。美国和欧洲价格，西德克萨斯原油价格和石油输出国组织原油价格。我们将开发一个回归模型，并尝试预测布伦特原油，美国价格给定的所有系列在我们的模型。

plt.plot(df['opec'], label="opec")
plt.plot(df['euro'], label="Euro")
plt.plot(df['br_crude'], label="Brent Crude Oil")
plt.plot(df['wti'], label="West Texas")plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.title('Opec Oil, Europe Brent Crude, Brent Crude Oil & West Texas Price over Time')
plt.show()

我们可以看到，所有的序列都遵循一个随机趋势，表现出更强的跨时变化和更大的下降；更重要的是，所有的系列似乎都有某种联系。这两个系列的水平看起来都不稳定。它们似乎有一个共同的趋势，表明它们可能是共同整合的。

IQR 绘制所有变量:

plt.title('IQR plot for all variables')
sns.boxplot(x=df)
plt.show()

我们可以在数据中看到一个小的异常值，可以忽略不计。

VAR 模型描述了平稳变量之间的动态相互关系。因此，时间序列分析的第一步应该是单位根检验，以确定序列是否平稳。

相关性检查:

pearsoncorr = df.corr(method='pearson')
sns.heatmap(pearsoncorr,
xticklabels=pearsoncorr.columns,
yticklabels=pearsoncorr.columns,
cmap='RdBu_r',
annot=True,
linewidth=0.5)
plt.show()

尽管相关性经常用于所有多变量金融时间序列用例，但是相关性可能非常不稳定。此外，ADF 或 Hurst 指数有助于我们从统计上确认该序列是否均值回复。然而，我们无法从这些测试中确定形成线性组合所需的套期保值比率，它们只会告诉我们，对于特定的情况，线性组合是否是稳定的。另一种安排是协整，这可能是两个金融序列之间联系的一个更稳健的措施。

增强型迪基-富勒(ADF):

ADF 检验表明，如果一个价格序列具有均值回复，那么下一个价格水平将与当前价格水平成比例。在数学上，ADF 是基于测试自回归时间序列样本中单位根的存在的思想。

def adfuller_test(series, signif=0.05, name='', verbose=False):
  r = adfuller(series, autolag='AIC') output = {'test_statistic':round(r[0], 4), 'pvalue':
round(r[1], 4), 'n_lags':round(r[2], 4), 'n_obs':r[3]}
  p_value = output['pvalue']
def adjust(val, length= 6): 
  return str(val).ljust(length)print(f'Augmented Dickey-Fuller Test on "{name}"', "\n   ", '-'*47)
print(f'Null Hypothesis: Data has unit root. Non-Stationary.')
print(f'Significance Level = {signif}')
print(f'Test Statistic = {output["test_statistic"]}')
print(f'No. Lags Chosen = {output["n_lags"]}')for key,val in r[4].items():
  print(f' Critical value {adjust(key)} = {round(val, 3)}') if p_value <= signif:
    print(f" => P-Value = {p_value}. Rejecting Null Hypothesis.")
    print(f" => Series is Stationary.")
  else:
    print(f" => P-Value = {p_value}. Weak evidence to reject the Null Hypothesis.")
    print(f" => Series is Non-Stationary.")

没有一个统计数据接近 5%的显著水平。因此，我们可以自信地说，级数在水平上是非平稳的。研究了它们的协整关系，在每种情况下(p 值> 0.05)，不能在任何合理的显著性水平上拒绝非平稳性的零假设。

协整检验:

协整检验过程中的步骤:

1. 使用单变量单位根检验，比如 ADF、PP 检验，分别检验每个成分序列的单位根。
2. 如果单位根不能被拒绝，那么下一步就是检验成分之间的协整，即检验αYt 是否为 I(0)。

如果我们发现时间序列是一个单位根，那么我们就进入协整过程。协整检验主要有三种方法:Johansen 方法、Engle-Granger 方法和 Phillips-Ouliaris 方法。我们将使用恩格尔-格兰杰检验。通过使用协整检验，我们将检验 wti 价格指数和 br_crude 价格指数之间是否存在协整的长期关系。数学上，我们将考虑协整关系:

Engle & Granger 的两步残差检验程序首先估计一个横截面回归，然后使用具有修正临界值的 ADF 分布检验该回归的残差。

这里记录了每个时间点的观测向量。在这种情况下，可能有一些构成平稳时间序列的向量的线性组合和其他非平稳的线性组合。

从上面的模型残差图来看，平均值为 0，这可能是一个趋势。

估计的协整向量非常接近[+1，1]，表明存在简单的无套利关系。让我们用一个常数和一个时间趋势(“CT”)重新运行测试。我们在 5% (-3.79)和 10% ( 3.50)水平的临界值的测试统计为-4.231<；p 值< 0.05，我们可以拒绝无协整的零假设。

残差明显为零，但显示了当前年份(2020 年)的结构性突变的证据。

同样，我们测试了其他可能的配对组合，以获得类似的输出。

虽然 Engle-Granger 方法很容易使用，但是它的一个主要缺点是它只能估计变量之间的一个协整关系。

这里，我们有 1 个以上的变量，可能有一个以上的线性独立的协整关系。

p 值的重要性:

现实生活中商业案例的 p 值一直是争论的焦点。p 值在发表的文章上可能看起来不错，但它从来没有像今天这样被使用过。更多详情请阅读这篇有趣的 [文章](http://The P value was never meant to be used the way it’s used today.) 。传统上， p- 值< 0.05 是显著性的标准。

OLS 回归:

现在，我们将使用 OLS 拟合时间序列之间的回归，并可视化残差图。如果残差是平稳的，我们可以有把握地假设，给定的序列确实是协整的。

residuals = results.resid
br_crude_fit = results.fittedvaluesfig = plt.figure(1)
ax1 = fig.add_subplot(111)
ax1.plot(br_crude_fit, label='Linear Prediction')
ax1.plot(df['br_crude'], label='Brent Crude-USA')
ax1.set_xlabel('Date')
ax1.legend(loc=0)
ax2 = plt.twinx()
ax2.set_ylabel('Residuals')
ax2.plot(residuals, label='Residuals')
ax2.legend(loc=0)
plt.grid(True)
plt.show()

残差的级别(红线)，它看起来比原始序列(蓝线)更像一个平稳序列。在这里，我们可以看到，实际和拟合线是难以区分的。

现在让我们对残差序列进行 ADF 检验。我们使用 DFGLS 测试，并指定 12 个滞后作为 DFGLS 回归的最大滞后阶数。

滞后的数量表明最佳滞后长度为 7。我们的测试统计(–7.059)< Critical Values at 1%, 5% & 10% and null hypothesis of no cointegration can be rejected.

breus ch-Godfrey(BD)测试:

BD 比 DW 更通用，允许我们测试高阶自相关。

残差的分布看起来像一个钟形，尽管有一些可能导致偏斜的异常值。

Jarque-Bera 正态性检验:

χ2 检验的双尾 p 值等于零意味着我们接受残差的零假设总体上是正态分布的，即峰度和偏度都是正态分布的。

向量自动回归(VAR):

VAR 方法使用 AR 模型对每个时间序列的下一步进行建模。模型的符号包括将 AR§模型的阶数指定为 VAR 函数{VAR§}的参数。

让我们先来看看这个系列的第一个不同点，然后再试试 ADF。通常，如果水平时间序列不是平稳的，第一个差异将是。数据差异或多或少是计算时间序列从一个值到下一个值的变化量。

df = df [['br_crude', 'wti', 'opec', 'euro']]
nobs = 5
df_train, df_test = df[0:-nobs], df[-nobs:]print(df_train.shape)
print(df_test.shape)df_tr_diff = df_train.diff()
df_tr_diff.dropna(inplace=True)
print(df_tr_diff.head())df_tr_diff.plot(figsize=(10,6), linewidth=5, fontsize=20)
plt.show()

考虑到时间因素，为时间序列问题创建一个训练和测试集是很棘手的。在这个阶段任何错误的举动都会打乱这个系列的格局。创建测试集时应该考虑日期和时间值。

现在，我们对数据帧的第一差重复上述所有步骤。

# ADF Test on each column
for name, column in df_tr_diff.iteritems():
adfuller_test(column, name=column.name)

这里，我们看到所有的级数在一阶差分后都是平稳的。

风险值模型:

var_model = smt.VAR(df_tr_diff)
res = var_model.select_order(maxlags=10)
print(res.summary())

#Fit VAR model
var_model = VAR(endog=df_tr_diff)
var_model_fit = var_model.fit(maxlags=10)#Print a summary of the model results
var_model_fit.summary()

我们在模型总结中真正想要关注的是 br_crude 的等式，其中 br_crude 根据自身的滞后值和其他变量的滞后价格来估计其他价格。t 统计值越高，我们就越有可能拒绝 H0，两个变量之间也越有可能存在相关性。p 值也用于拒绝零假设。p 值越低，反对零假设的证据越强。

估计 VAR 的方程式并不需要很强的假设；然而，计算脉冲响应函数(IRF)或方差分解确实需要识别限制。一个典型的限制是假设一对变量之间的动态关系，例如，x 只滞后地影响 y，或者 x 在长期内不影响 y

让我们使用脉冲响应图，在不同的水平上，将一个变量的变化对其他变量的影响可视化。

# Impulse Response Analysis
irf = var_model_fit.irf(20)
irf.plot()
plt.show()

给定参数估计和 Engle Granger 测试结果，这里建立了序列之间的联系。可以看出，对冲击的响应在这里被捕获，并且它们几乎在第 10 个滞后之后消失。

预测误差差异分解(FEVD):

在 FEVD 曲线图的第一行中显示了对 br_crude 价格的冲击对其它价格和它们自身的其它价值的影响。在 br _ 原油价格(深黑色柱)的情况下，归因于自身冲击的误差百分比为 100%，因为 wti 系列解释了约 40%的回报变化；欧佩克约为 60%，欧元系列解释了约 70%的变化。

注:上述测试程序的概念取自 Brooks(2019)；参见参考部分。

预测:

现在模型已经建立，让我们生成预测并与测试/验证集中的实际数据进行比较。

# Get the lag order
lag_order = var_model_fit.k_ar
print(lag_order)# Input data for forecasting
input_data = df_tr_diff.values[-lag_order:]
print(input_data)# forecasting
pred = var_model_fit.forecast(y=input_data, steps=nobs)
pred = (pd.DataFrame(pred, index=df_test.index, columns= df.columns + '_pred'))
print(pred)

反转变换:

def invert_transformation(df_tr_diff, pred):
  forecast = pred.copy()
  columns = df_train.columnsfor col in columns:
  forecast[str(col)+'_pred'] = df_train[col].iloc[-1] +  
                               forecast[str(col)+'_pred'].cumsum()
  return forecastoutput = invert_transformation(df_tr_diff, pred)
output.loc[:, ['br_crude_pred']]

预测评估:

# forecast bias
forecast_errors = [combine['br_crude'][i]- combine['br_crude_pred'][i] for i in range(len(combine['br_crude']))]
bias = sum(forecast_errors) * 1.0/len(combine['br_crude'])
print('Bias: %f' % bias)# MAE
mae = mean_absolute_error(combine['br_crude'], combine['br_crude_pred'])
print('MAE: %f' % mae)# MSE & RMSE
mse = mean_squared_error(combine['br_crude'], combine['br_crude_pred'])
print('MSE: %f' % mse)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)

实际与预测:

combine = pd.concat([df_test['br_crude'], output['br_crude_pred']], axis=1)
combine['accuracy'] = round(combine.apply(lambda row: row.br_crude_pred /
row.br_crude *100, axis = 1),2)combine['accuracy'] = pd.Series(["{0:.2f}%".format(val) for val in combine['accuracy']],index = combine.index)
combine = combine.assign(day_of_week = lambda x: x.index.day_name())
combine = combine.round(decimals=2)
combine = combine.reset_index()
combine

fig = go.Figure()
fig.add_trace(go.Scatter(x=combine['Date'], y=combine['br_crude'], name="Actual Crude price"))
fig.add_trace(go.Scatter(x=combine['Date'],y=combine['br_crude_pred'],name="Predicted crude price"))
fig.update_layout(
title="Actual vs Predicted Brent Crude Price",
xaxis_title="Date", yaxis_title="Price", font=dict(family="Courier New, monospace",size=18,color="#7f7f7f"))
fig.update_layout(autosize=False,width=1000,height=400,)
fig.update_layout(legend_orientation="h")
fig.show()

结论:

一旦我们有了数据的模型，分析如何评估它的质量是很重要的。第一种选择是使用残差。残差基本上是预测值和实际值之间的平方差。这可能是一个简单的选择来判断我们的模型的准确性。

残差与我们用来创建模型的数据量成比例。因此，不建议使用所有残差的和。一种选择是将残差除以瞬时数量乘以因变量的数量。

我这里可以到达。

参考文献:

1. 汉密尔顿法学博士，1994 年。时间序列分析。新泽西，普林斯顿大学出版社。
2. 恩格尔，R.F .，格兰杰，C.W.J，1987。协整与误差修正:表示、估计与检验。计量经济学。55, 251–276.
3. 布鲁克斯，C. (2019 年)。金融计量经济学导论。剑桥大学出版社。