import numpy as npimport matplotlib.pyplot as plt% matplotlib inline
本章给出的例子与Cameron Davidson Pilon在2015年出版的Bayesian Methods for Hackers: Probabilistic Programming and Bayesian Inference一书中所给的例子类似。但是，我们的例子被重新编写以更加适配金融背景，书中的数学观点也可以从代码中更直不雅观的表示。
假设你有一种证券，它要么值1美元，要么一文不值。收益取决于两步：（1）证券的收益随机，其概率为50%。在收益随机的情形下，50%可能得到1美元收益，50%收益为0；（2）证券拥有真正的收益，其概率也为50%。在这种情形下，得到1美元收益的概率
，称为真正收益概率（True Payoff Probability，TPP）。
这个收益方案如图10.1所示。
图10.1 收益方案
你更感兴趣的是找出真正的收益率是多少，由于真正的收益率将会见告你采纳什么交易策略。在本例中，你购买100股证券，在这100股证券中，有54股让你赚了1美元。
实际的TPP是多少呢？在本例中，可以通过给出解析解来剖析最可能的TPP。但是我们将利用一种适用于更繁芜问题的打算方法。
在第10.1.1节中，我们将证券收益过程进行仿照。
10.1.1 扁平先验变量
表示TPP。我们随机采样100个真值；如果你在真正收益情形下得到了1美元，则该真值为1，否则为0。我们在图10.1的受益方案中，在“开始”和“随机收益率”两个过程也进行随机抽样。只管有的时候并不是全都须要，但对所有试验一口气全部采样的随机结果在打算效率上更高效。
末了，我们将所有收益相加，然后除以证券的数量，终极我们求的是本次仿真试验中的收益率。
下面的代码段是仿真运行函数。然而，你须要确保能够理解打算过程是如何与图10.1的方案对应的，这是非常主要的事情。[1]
def run_sim(x): truth = np.random.uniform(size=100) < x first_random = np.random.randint(2,size=100) second_random = np.random.randint(2,size=100) res = np.sum(first_randomtruth +(1-first_random)second_random)/100 return res
接下来，我们想试验一些可能的TPP。因此，本例将对一个候选TPP进行抽样，并利用此概率进行仿照。如果仿真输出的收益与我们在现实生活中不雅观察到的收益相同，那么这个候选TPP便是一个真实概率。
下面的采样函数将返回真实概率。当候选概率与输入不等的时候，返回None：
def sample(data = 0.54): x = np.random.uniform() if run_sim(x) == data: return x
由于须要对可能的TPP进行抽样，因此我们很自然地希望能加快这一过程。为此，我们利用JobLib的库，采取并行实行的办法来加速上述打算过程。
提示：
JobLib在Kaggle内核上已经预安装。
我们须要导入Parallel类和delayed方法。Parallel类将有助于并行地实行这些循环，delayed方法则有助于在并行循环内按顺序实行相应函数。通过以下代码来导入它们：
from JobLib import Parallel, delayed
上述Parallel类和delayed方法的详细细节与本例无关。Parallel(n_jobs=-1)方法使任务的并行实行数量与机器上的CPU数量相同。例如，delayed(sample)() for i in range(100000) 代码中sample方法将实行100000次。[2]
我们得到Python列表t，并将其转换为一个NumPy数组。正如你在以下代码段中看到的，数组中大约有98%的值为None。这意味着在采样器测试的
值中，有98%与我们的数据不匹配。
t = Parallel(n_jobs=-1)(delayed(sample)() for i in range(100000))t = np.array(t,dtype=float)share = np.sum(np.isnan(t))/len(t)100print(f'{share:.2f}% are throwaways')98.01% are throwaways
因此，我们现在将丢弃所有None值，留下
值：
t_flat = t[~np.isnan(t)]plt.hist(t_flat, bins=30,density=True)plt.title('Distribution of possible TPPs')plt.xlim(0,1);
运行上述代码后，我们将得到如图10.2所示的输出。
正如你所看到的，图10.2展示了可能的TPP分布情形。这张图显示，TPP最可能的范围是50%～60%。虽然其他值也是有可能的，但它们的可能性更小。
你刚才看到的是贝叶斯方法的优点之一。所有的估计以分布形式展示，我们可以打算置信区间（在贝叶斯术语中称为“可信区间”）。
图10.2 朴素采样器创造的可能真实收益概率的分布
这让我们能够对事情的确信度更加精确，也让模型中其他参数值的可能性更加精确。将它与我们感兴趣的金融领域联系起来。在金融运用中，数以百万美元的金钱都按照模型输出进行交易，因此能够量化这种不愿定性就是非常大的上风。
10.1.2 < 50%先验此时，你可以将结果提交给证券交易领域的一名专家，他看了你的剖析，摇了摇头说：“TPP不能超过0.5。”他阐明道，“从根本业务来看，这在客不雅观上不可能超过50%。”
那么，你将如何把这个事实领悟到仿真剖析中呢? 最直接的办理方案是仅试验从0到0.5的候选TPP。你所须要做的便是限定候选
的采样空间。可以通过运行以下代码来实现：
def sample(data = 0.54): x = np.random.uniform(low=0,high=0.5) if run_sim(x) == data: return x
现在，你可以像刚才那样进行仿真了：
t = Parallel(n_jobs=-1)(delayed(sample)() for i in range(100000))t = np.array(t,dtype=float)# Optionalshare = np.sum(np.isnan(t))/len(t)100print(f'{share:.2f}% are throwaways')99.10% are throwawayst_cut = t[~np.isnan(t)]plt.hist(t_cut, bins=15,density=True)plt.title('Distribution of possible TPPs')plt.xlim(0,1);
像前面那样，上面代码将输出图10.3所示的结果。
图10.3
从0～0.5区间下潜在TPP分布
10.1.3 先验与后验显然，你选择的值会影响仿真剖析的结果，这也反响出你对
可能的值的认识。
在第一轮中，在你看到任何数据之前，认为所有介于0和100%之间的TPP都是等概率的，这被称为“扁平先验”，由于所有值的分布都是相同的，因此是扁平的。在第二轮中，你认为TPP必须低于50%。
在看到数据之前能表示你对
认识的分布称为“先验分布
”，简称“先验”。我们从仿真中得到的可能
的分布（也便是看到数据
后）称为“后验分布
”，简称“后验”。
以下代码给出了第一轮和第二轮的先验和后验采样情形。首先是绘制扁平后验的结果，运行如下代码：
flat_prior = np.random.uniform(size=1000000)plt.hist(flat_prior,bins=10,density=True, label='Prior')plt.hist(t_flat, bins=30,density=True, label='Posterior')plt.title('Distribution of $x$ with no assumptions')plt.legend()plt.xlim(0,1);
上述代码天生图10.4。
图10.4 扁平先验的采样结果
接下来是对小于50%先验的采样输出。运行如下代码：
cut_prior = np.random.uniform(low=0,high=0.5,size=1000000)plt.hist(cut_prior,bins=10,density=True, label='Prior')plt.hist(t_cut, bins=15,density=True, label='Posterior')plt.title('Distribution of $x$ assuming TPP <50%')plt.legend()plt.xlim(0,1);
上述代码天生图10.5。只管利用了相同的采样器，但你可以看到结果是非常不同的。
图10.5 对小于50%先验的采样结果
你把稳到什么奇怪的事情了吗?第二轮的后验值大致即是第一轮的后验值，但第二轮后验值在0.5处被截断了。这是由于第二轮试验中当
大于0.5的时候，它的先验值为0，对付其他值为1。
由于我们只记录与数据相匹配的仿真结果，因此直方图中显示的仿真结果反响了在给定TPP和
的情形下，运行仿真试验并输出不雅观察数据
的概率
。
从仿真试验得到的后验概率
即是在给定某个TPP情形下进行试验并不雅观察到数据的概率
乘以概率
。
数学上，上述关系表示为：
当数据可以通过诸如面对面会议等办法自然地获取时，那么我们可能须要解释数据网络方法中存在的偏见。大多数情形下，我们不必担心这个问题，可以大略地忽略它。但有些时候，这种丈量能够放大某些结果。
为了缓解这个问题，我们须要除以数据分布
将其作为后验公式的末了补充，并得到如下公式：
如你所见，这便是贝叶斯公式！
当我们做仿真试验时，是从后验中采样。为什么不能直接用贝叶斯公式来打算后验呢？答案很大略，打算
须要我们对TPP进行积分。这个打算是非常棘手的。我们的仿真方法是一种大略方便的替代办理方案。
提示：
第一轮先验（即所有的TPP都是等可能的）称为扁平先验，由于我们对TPP值的分布不做任何假设。在这种情形下，贝叶斯后验即是最大似然估计。
10.1.4 马尔可夫链蒙特卡罗算法在第10.1.3节中，我们通过从先验中随机采样，然后试验采样值来近似估计后验分布。这种随机试验的方法在模型只有一个参数（比如TPP）时效果会很好。然而，随着模型繁芜度的提高和参数的增加，随机搜索方法将变得更慢。终极，将可能有太多的参数组合没有机会天生数据。因此，我们须要频繁地辅导搜索操作并利用更高的后验概率来采样参数。
这种有指向性但仍旧随机的采样方法被称为马尔可夫链蒙特卡罗算法。蒙特卡罗表示随机和仿真，马尔可夫链表示我们在参数空间上以某些概率进行搜索。
在这里谈论的特定算法中，我们将以一定的概率采取一个不同的参数值，这个概率是该参数值的后验概率的比率。这里，我们将考虑打算参数值的后验概率。由于概率不能大于1，我们将把比率限定在1以内。但这只是一个数学上的限定，对算法来说并无太大影响。
图10.6展示了马尔可夫链蒙特卡罗算法的基本事情事理。
图10.6 马尔可夫链蒙特卡罗算法
图10.6显示我们处于“随机行走”中，或多或少地随机遍历不同的参数值。然而，实际上我们并不是完备“随机移动”，而是更方向于具有更高后验概率的参数值。
为了实行这个算法，我们须要做4件事。
1．从当前参数值
处提出一个新的参数值
。
2．利用贝叶斯法则来估计参数
的后验概率
。
3．打算从
移动到新参数值
的概率
。（把稳概率必须小于1。）
4．以概率
转移到新的参数值。
接下来，我们逐步构建这个算法的模块：
# REPETITION FROM FIRST SECTIONdef run_sim(x): truth = np.random.uniform(size=100) < x first_random = np.random.randint(2,size=100) second_random = np.random.randint(2,size=100) res = np.sum(first_randomtruth +(1-first_random)second_random)/100 return res# REPETITION FROM FIRST SECTIONdef sample(x,data = 0.54): if run_sim(x) == data: return x
首先，我们须要设置一个新的
。由于我们不想盲目地随机搜索，而是更精确地随机游走，因此设置
须要依赖于
的先前值。在本例中，我们将从均值为
、标准差为0.1的正态分布中采样
。
只要知足
与
是干系的条件，从其他分布或具有不同标准差的正态分布中取样也是可以的：
def propose(x): return np.random.randn() 0.1 + x
在本章第一部分内容中，我们通过从先验中采样然后通过运行仿真程序的方法来实现直接从后验中采样。现在，我们通过马尔可夫链蒙特卡罗算法来采样，这样就不再直接从后验采样。因此，可以利用贝叶斯法则来打算后验概率。
记住常日不除以
，由于我们不假设有偏丈量。贝叶斯法则可以简化为
，个中
为后验概率，
为先验概率，
为似然值。因此，为了估计参数
的似然值，我们利用该参数进行了多次仿真试验。
似然值是仿真试验匹配真实数据的比例：
def likelihood(x): t = Parallel(n_jobs=-1)(delayed(sample)(x) for i inrange(10000)) t = np.array(t,dtype=float) return (1 - np.sum(np.isnan(t))/len(t))
我们首先再次利用扁平先验打算，也便是每个TPP都是等概率的：
def prior(x): return 1 #Flat prior
参数
的后验概率为似然值乘以先验概率：
def posterior(x): return likelihood(x) prior(x)
现在，我们准备将上面所有内容整合到Metropolis-Hastings MCMC算法中。
首先，我们为
设置初始值。为了让算法快速找到可能的值，明智的做法是将
初始化为最大似然值或我们认为可能的某个估计值。如果要打算这个初始值的后验概率，运行以下代码即可：
x = 0.5pi_x = posterior(x)
同样，须要记录试验中采样的所有值。为了展示，我们也跟踪记录了后验概率。运行以下代码：
trace = [x]pi_trace = [pi_x]
现在，我们进入主循环。在此之前，我们须要记住算法包含4个步骤：
1．设置新的候选
;
2．打算后验概率
;
3．打算接管概率：
4．以概率
来设置
为
。
for i in range(1000): #Main Loop x_cand = propose(x) pi_x_cand = posterior(x_cand) alpha = np.min([1,pi_x_cand/(pi_x + 0.00001)]) # Save division u = np.random.uniform() (x, pi_x) = (x_cand,pi_x_cand) if u<alpha else (x,pi_x) trace.append(x) pi_trace.append(pi_x) if i % 10 == 0: print(f'Epoch {i}, X = {x:.2f}, pi = {pi_x:.2f}')Epoch 0, X = 0.50, pi = 0.00Epoch 10, X = 0.46, pi = 0.04...Epoch 990, X = 0.50, pi = 0.06g
在运行这个算法多少轮后，我们终极得到一个潜在作弊者收益份额的分布。跟之前所做的一样，运行以下代码使结果可视化：
plt.hist(trace,bins=30)plt.title('Metropolis Hastings Outcome')plt.xlim(0,1);
结果如图10.7所示。
图10.7 Metropolis Hastings采样器的输出
通过查看trace随韶光变革的轨迹，trace数据显示了算法是如何以高度可能的值为中央进行随机移动的：
plt.plot(trace)plt.title('MH Trace');
我们将以图脸色势得到输出结果如图10.8所示，该结果向我们展示了Metropolis Hasings采样器的运行轨迹。
图10.8 Metropolis Hastings采样器的轨迹
为了更好地理解，我们绘制了后验概率随所试验值变革的情形：
plt.scatter(x=trace,y=pi_trace)plt.xlabel('Proposed X')plt.ylabel('Posterior Probability')plt.title('X vs Pi');
成功实行代码后，我们将得到图10.9所示的输出图像。
图10.9 试验值和后验概率
10.1.5 Metropolis-Hastings MCMC为了展示PyMC3的强大功能和灵巧性，我们将它用于一个经典的计量经济学问题，但重点要从贝叶斯角度对其进行阐述。
提示：
本例直接改编自PyMC3文档中的示例，而PyMC3文档中的例子改编自霍夫曼论文“No-U-Turn Sampler”。
股票价格和其他金融资产价格都颠簸，日收益的方差被称为颠簸率。颠簸率是一种常用的风险评估方法，因此对其进行准确的度量是非常主要的。
大略的办理方案是打算收益的事后风险度量，即方差。但是，表达实际收益颠簸率的不愿定性是有益的。与我们之前看过的收益例子类似，这里存在一个实际值的分布，从这个分布中可以得到已实现的值。由于存在可能的颠簸率值的分布，个中不雅观察到的颠簸率可能便是一个已实现的样本，这也称为“随机颠簸率”。
在本例中，我们建立标准普尔500指数的随机颠簸率模型。为此，我们首先加载数据，你可以直接从互联网下载数据，也可以在Kaggle平台上找到数据。
运行以下代码加载数据：
df = pd.read_csv('../input/S&P.csv')df['Date'] = pd.to_datetime(df['Date'])
在这个例子中，我们对收盘价感兴趣，因此须要从数据集中提取收盘价。数据集先显示新的数据。因此须要将其反转，通过以下代码实现：
close = pd.Series(df.Close.values,index=pd.DatetimeIndex(df.Date))close = close[::-1]
以下代码绘制了收盘价。
close.plot(title='S&P 500 From Inception');SP500
我们将得到图10.10所示的图表。
图10.10 从股市开市以来到2018年S&P500收盘价
数据集包含了标准普尔自开市以来的数据，这对我们来说有点太多了。以是在本案例中，我们将其在1990年截断。通过运行下列代码来指定这个韶光：
close = close['1990-01-01':]
由于我们对收益感兴趣，因此须要打算价差。我们利用np.diff来获取逐日价差。为了便于画图，我们将逐日价差结果用pandas Series工具来存储：
returns = pd.Series(np.diff(close.values),index=close.index[1:])returns.plot();
输出结果如图10.11所示。
图10.11 从1990年到2018年的S&P500的收益
PyMC3包含一些处理韶光序列所需的分外分布，如随机游走。当我们想要对股票价格进行建模时，利用PyMC3是非常精确的方法。
首先，导入PyMC3及其时间序列工具random walk类：
import pymc3 as pmfrom pymc3.distributions.timeseries import GaussianRandomWalk
然后，建立模型，可以运行以下代码来实现建模：
with pm.Model() as model: step_size = pm.Exponential('sigma', 50.) #1 s = GaussianRandomWalk('s', sd=step_size, #2 shape=len(returns)) nu = pm.Exponential('nu', .1) #3 r = pm.StudentT('r', nu=nu, #4 lam=pm.math.exp(-2s), observed=returns.values)
现在让我们回顾下建模须要实行的代码命令。如你所见，它有4个关键要素组成。
1．颠簸率
被建模为随机游动过程，该随机游走的步长为step_size。步长的先验分布是
= 50的指数分布。（再次解释，理解所利用的每个分布的细节知识对付演示和试验来说并不是必需的。）
2．然后对随机颠簸率本身进行建模。由于步长本身便是一个随机变量，请把稳我们是如何将步长融入随机颠簸率模型中的。随机游走的长度应与不雅观察到的收益值相同。
3．我们建立了知足nu个自由度的StudentT分布的实际股票收益模型。nu的先验也是指数分布。
4．末了，我们对实际收益进行建模。将实际收益建模成比例因子为lam的StudentT分布，个中lam是由随机颠簸模型天生的。为了使模型基于不雅观察到的数据进行构建，我们将不雅观察的收益值作为参数传入模型中。
PyMC3的标准采样器不是Metropolis Hastings，而是No-U-Turn采样（No-U-Turn Sampler，NUTS）。如果我们不指定采样器而仅调用sample，则PyMC3将默认为NUTS。
为了使采样平稳运行，须要指定较大量的tune采样。采样器将会从这些tune采样中抽取样本，以便找到一个好的出发点。类似于之前的burned采样，这些tune采样也不会成为后验中的一部分。
我们须要设置较高的target_accept值，来见告NUTS在接管值时要宽容。可以运行以下代码来实现：
with model: trace = pm.sample(tune=2000,nuts_kwargs=dict(target_accept=.9))
PyMC3有一个很好用的工具，它可用于可视化采样结果。我们对颠簸率随机游动的标准偏差
以及实际收益所符合的StudentT分布的自由度比较感兴趣。
在并走运行两段代码时，你可以看到两个不同的输出分布。如果我们将采样器运行更长的韶光，那么这两个结果将会收敛。可以通过对它们进行均匀来得到更好的估计，这便是PyMC3针对预测所采取的方法。现在让我们考试测验利用以下代码：
pm.traceplot(trace, varnames=['sigma', 'nu']);TracePlot
图10.12显示了该代码的实行结果。
图10.12 PyMC3采样器的结果概览（左侧两个采样器天生的分布，右侧是它们的轨迹）
在末了一步，我们展示了随机颠簸率是如何随韶光变革的。你可以看到它是如何与经济颠簸周期（例如2008年金融危急）完美契合的。你还可以看到，在某些期间中，模型对颠簸率大体是确定的：
plt.plot(returns.values)plt.plot(np.exp(trace[s].T), 'r', alpha=.03);plt.xlabel('time')plt.ylabel('returns')plt.legend(['S&P500', 'Stochastic Vol.']);
该代码的输出如图10.13所示。
图10.13 1990年至2018年的推断随机颠簸率
有大量的运用可以利用此类相对小规模的贝叶斯模型来进行建模，其紧张优点是模型易于阐明并且可以很好地表达不愿定性。由于模型很清楚地表述了需求，因此概率编程与数据科学中的“storytelling”方法很同等。
在第10.1.6节中，我们将从浅层概率编程转向深层概率编程。
10.1.6 从概率编程到深度概率编程到目前为止，我们开拓的贝叶斯模型都很浅。因此，我们要考虑是否可以将深度网络的预测能力与贝叶斯模型的上风相结合，这是一个生动的研究领域。
深度网络具有大量参数，这使得在参数空间中搜索成难堪题。在传统的监督深度学习中，我们利用反向传播来办理这个问题，反向传播也可以用于贝叶斯模型。但是，它不是进行贝叶斯深度学习的唯一方法，乃至不一定是最佳方法。
总体来说，有4种方法可以用于贝叶斯深度学习。
利用自动微分变分推理（Automatic Differentiation Variational Inference，AVI）。这意味着用勾引模型来近似后验分布，然后利用梯度低落来优化模型参数。PyMC3利用AVI优化器可以实现这一点。请参阅Alp Kucukelbir等人于2016年揭橥的论文“Automatic Differentiation Variational Inference”。或者，你可以利用Pyro，它实现了快速的、GPU优化的AVI。只管在这里不适宜花大量篇幅先容PyMC3，但还要推举PyMC3文档，它是一个很好的教程。假设后验值呈正态分布，然后利用标准神经网络库（例如Keras）并学习每个参数的均值和标准差。还记得我们在利用变分自编码器时如何从参数化正态分布中采样值吗？我们可以对每一层都这样做。与AVI比较，这种方法演习速率更快，占用的打算资源和内存更少，但灵巧性较差，其参数是非贝叶斯神经网络的两倍。利用dropout技巧。在利用韶光序列时，我们在测试韶光打开了dropout功能，并多次运行推理以得到置信区间。这是一种非常随意马虎实现的贝叶斯学习，没有比常规神经网络更多的参数。但是，它在推理时速率较慢，并且也不具有AVI的所有灵巧性。挑选并稠浊。为了演习神经网络，我们须要一个可以从AVI得到的梯度旗子暗记。我们可以以常规办法演习神经网络的socket（有时我们称之为“特色提取器”），并以贝叶斯办法演习网络的头部。这样，我们可以得到不愿定性估计，而不必承担贝叶斯方法的全部本钱。本文摘自：《金融中的机器学习》
机器学习是设计与运用算法的科学，可从数据中进行学习和预测，其运用已经非常普遍。金融领域集中了大量的交易数据，为人工智能技能的利用奠定了良好的数据根本。本书面向金融领域的读者，先容了机器学习技能的事理与实践。
本书包括10章，先容了神经网络算法、构造化数据的处理、打算机视觉处理技能、韶光序列剖析、自然措辞处理、天生模型的运用、强化学习技能、数据建模与调试、贝叶斯推理和概率编程等内容。
本书由资深金融从业者编写，领悟了其在金融项目中关于机器学习的实践履历，适宜金融领域的数据科学家、数据剖析师、金融科技公司的技能研发职员以及对金融领域的机器学习技能感兴趣的读者阅读。