手把手：基于概率编程Pyro的金融预测，让正则化结果更有趣！

大数据文摘作品

编译：修竹、笪洁琼、夏雅薇

作者用了一种新奇的方法来训练神经网络。更新权重的分布而不是顺序更新静态权重，得到了更有趣和可靠的结果。贝叶斯方法给了我们一个机会，使得我们可以不手动添加正则项的情况下对神经网络进行正则化，理解模型的不确定性，并尽可能使用更少的数据得到更好的结果。

Hi！又见面啦。去年我推出了几篇基于神经网络的金融预测教程，我认为有些结果还是蛮有趣的，值得应用在实际交易中。

如果你读过那些教程，你一定会注意到，当你试图在“随机”数据上用一些机器学习模型并且希望找到隐藏模式时，你其实正逐渐对训练集进行过拟合。

我们使用不同的正则化方法和补充数据来解决这个问题，但是这非常耗时间并且有点盲目搜索了。

今天我想介绍一种稍微不同的方法来用于相同的算法。从概率角度讲，我们可以从数据本身学习正则化方法，在我们预测中估计准确性，使用更少的数据来训练并且在模型中加入概率依赖。

我不会深入到贝叶斯模型或变分推理的技术或者数学细节上，我将给出一些概述，同时也会更加关注如何应用。像往常一样，你可以在下面的链接内查看代码。

代码链接：

https://github.com/Rachnog/Deep-Trading/tree/master/bayesian

为了更深入了解概率编程、贝叶斯模型及其应用，我推荐以下资源给大家：

• 模式识别和机器学习：

  http://www.springer.com/us/book/9780387310732

• 为黑客设计的贝叶斯方法：

  https://www.amazon.com/Bayesian-Methods-Hackers-Probabilistic-Addison-Wesley/dp/0133902838

同时推荐以下python库：

• PyMC3：

  https://github.com/pymc-devs/pymc3

• Edward：

  http://edwardlib.org/

• Pyro：

  http://pyro.ai/

概率编程

这个概率性的东西是什么，而且我们为什么要称之为编程呢？首先，我们先回忆一下“正常”的神经网络以及我们能从中获得什么。

我们有参数（权重），这些参数以矩阵表示，输出通常是一些标量值或者向量（例如用于分类时）。比如说，在用SGD训练模型之后，我们有了这些固定矩阵和网络在相同的输入样本上输出相同的向量。完全正确！

但是如果我们认为这些参数和输出都是互相依赖的分布呢？

神经网络中的每个权重都是来自某个分布的样本，输出也一样，每个输入来自整个网络的样本，同时这个网络依赖参数的样本。它给予了我们什么？

我们从最基础的开始讲。

如果我们把网络看成一组彼此依赖的分布，首先定义联合概率分布为p(y, z|x)，输出为y，还有一些依赖输入x 的模型“内部”、隐藏参数z（和普通神经网络一样）。

而我们想要找到一种神经网络分布，我们可以对y ~ p(y|x)采样然后把分布作为输出（该分布的样本期望值通常就是输出，标准差用来评估不确定性，如果分布模型的尾部越大，我们对于输出越没有信心）。

这个设置或多或少已经很清楚了，我们只需要记住，现在所有的参数，不管是模型的输入还是输出，都是分布。我们需要的训练是找到这些分布的参数以便在实际任务中获得更高的准确率。

必须要提到的是，参数分布的形状是我们自己设置的（例如，所有的初始权重都是w ~ Normal(0, 1)，然后我们将学习正确的均值和方差）。

初始分布称之为先验分布，使用过训练数据拟合参数的分布叫做后验分布。后者用于取样和获得输出数据。

模型的拟合效果怎么样呢？一般的框架叫做变分推理。我们不去深入了解细节，在这里我们需要寻找的模型是可以最大化似然函数log p_w(z|x)的, w 是模型的参数（分布参数），z是隐藏变量（隐藏神经元输出，从参数为 w 的分布中取样得到的），x是输入数据样本。这就是我们的模型。

在Pyro库中我们引入了一个实例作为这个模型的指导，指导中包括一些对所有隐藏变量q_ф(z)的分布，其中 ф叫做变分参数。这个分布必须近似于拟合数据最好的模型参数的“真实”分布。

训练的目的是最小化一个指导中关于输入数据和样本[log(p_w(z|x)) — log(q_ф(z))] 的期望。我们这里不讨论训练过程的细节，因为这里面包含好几门大学课程，现在我们就做黑盒优化好了。

哦对了，为什么是编程呢？因为我们通常将这种概率模型（比如神经网络）描述为从一个变量到另一个变量的有向图，这样我们就可以直接表示变量的依赖性：

最初这种概率编程语言被用来定义这些模型并对其进行推断。

为什么用概率编程？

你可以将它认为是一种附加的隐藏变量，从数据中学习模型，而不是采用在模型中注入dropout或L1正则化的方法。

考虑到所有权重都是分布，你可以从中进行N次抽样然后得到输出的分布，通过标准差可以估算你的模型对于结果的准确性。

而且还有个不错的赠礼就是我们只需要用更少的数据来训练模型，并且我们可以在变量间灵活的增加不同的依赖关系。

为什么不用概率编程呢？

我对于使用贝叶斯模型没有太多经验，但就我从Pyro和PyMC3学习中可以知道，训练过程耗时很长而且很难定义准确的先验分布。此外，处理分布的多个样本会导致误解和歧义。

数据展示

我获取了以太坊每日价格。这里面包括典型OHLCV（高开低走）元组以及每天关于以太坊推特的数量。

以太坊价格来源：

https://bitinfocharts.com/

我们将使用7天的价格、交易量和推特数量变化的百分比来预测下一天变化的百分比。

价格、推特数量和交易量变化

上图所示是数据的样本——蓝色表示价格变化，黄色表示推特数量变化，绿色表示交易量变化。这些值（0.1-0.2）之间有一些正相关，所以我们希望利用一些数据训练模型。

贝叶斯线性回归

首先我想了解简单线性回归在我们任务中的表现。

Pyro官方教程：

http://pyro.ai/examples/bayesian_regression.html

我们在PyTorch中定义了我们模型（详细解释请看官方教程）：

class RegressionModel(nn.Module):
 def __init__(self, p):
 super(RegressionModel, self).__init__()
 self.linear = nn.Linear(p, 1)
def forward(self, x):
 # x * w + b
 return self.linear(x)

这只是我们曾经用过的一个简单确定性模型，但是这也是在Pyro中定义概率的方法。

def model(data):
 # Create unit normal priors over the parameters
 mu = Variable(torch.zeros(1, p)).type_as(data)
 sigma = Variable(torch.ones(1, p)).type_as(data)
 bias_mu = Variable(torch.zeros(1)).type_as(data)
 bias_sigma = Variable(torch.ones(1)).type_as(data)
 w_prior, b_prior = Normal(mu, sigma), Normal(bias_mu, bias_sigma)
 priors = {'linear.weight': w_prior, 'linear.bias': b_prior}
 lifted_module = pyro.random_module("module", regression_model, priors)
 lifted_reg_model = lifted_module()
with pyro.iarange("map", N, subsample=data):
 x_data = data[:, :-1]
 y_data = data[:, -1]
 # run the regressor forward conditioned on inputs
 prediction_mean = lifted_reg_model(x_data).squeeze()
 pyro.sample("obs",
 Normal(prediction_mean, Variable(torch.ones(data.size(0))).type_as(data)),
 obs=y_data.squeeze())

上述代码中我们为参数W和b设置了一般线性回归模型分布，均为 ~Normal(0, 1)。我们称之为先验，创建了Pyro的随机函数（在我们例子中，是PyTorch的回归模型），将先验概率加到({‘linear.weight’: w_prior, ‘linear.bias’: b_prior})并且基于输入数据x从模型p(y|x)中采样。

这个模型的指导如下所示：

def guide(data):
 w_mu = Variable(torch.randn(1, p).type_as(data.data), requires_grad=True)
 w_log_sig = Variable(0.1 * torch.ones(1, p).type_as(data.data), requires_grad=True)
 b_mu = Variable(torch.randn(1).type_as(data.data), requires_grad=True)
 b_log_sig = Variable(0.1 * torch.ones(1).type_as(data.data), requires_grad=True)
 mw_param = pyro.param("guide_mean_weight", w_mu)
 sw_param = softplus(pyro.param("guide_log_sigma_weight", w_log_sig))
 mb_param = pyro.param("guide_mean_bias", b_mu)
 sb_param = softplus(pyro.param("guide_log_sigma_bias", b_log_sig))
 w_dist = Normal(mw_param, sw_param)
 b_dist = Normal(mb_param, sb_param)
 dists = {'linear.weight': w_dist, 'linear.bias': b_dist}
 lifted_module = pyro.random_module("module", regression_model, dists)
return lifted_module()

这里我们定义了我们想要“训练”的分布的变分分布。就像你看到的，我们为W和b定义了相同形状的分布，但是尽量使他们更接近实际（只要我们能想到的）。在这个例子里，我选择让这个分布形状更窄一些。

（~Normal(0, 0.1)）

我们以这种方式训练模型：

for j in range(3000):
 epoch_loss = 0.0
 perm = torch.randperm(N)
 # shuffle data
 data = data[perm]
 # get indices of each batch
 all_batches = get_batch_indices(N, 64)
 for ix, batch_start in enumerate(all_batches[:-1]):
 batch_end = all_batches[ix + 1]
 batch_data = data[batch_start: batch_end]
 epoch_loss += svi.step(batch_data)

之后我们想从模型中取样y。重复取样100次然后计算每一次取样预测的均值和标准差（标准差越大，我们对预测准确的信心越低）。

preds = []
for i in range(100):
 sampled_reg_model = guide(X_test)
 pred = sampled_reg_model(X_test).data.numpy().flatten()
preds.append(pred)

我们应该记得，金融预测中MSE，MAE或者MAPE等经典指标可能会让人很困惑——相对较小的错误率并不意味着你的模型运行良好，所以在样本外的数据上查看性能是非常重要的，这也是我们要做的：

30天的贝叶斯模型预测

如上图所示，结果并不是很好，但是最后一条的预测形状很好，这给了我们一点信心。让我们继续吧！

普通神经网络

在这个非常简单的模型之后，我们想试着尝试些更有趣的东西，比如神经网络。首先让我们先了解一个简单的MLP，只有一个隐藏层，包括含有25个神经元以及线性激活模型。

def get_model(input_size):
 main_input = Input(shape=(input_size, ), name='main_input')
 x = Dense(25, activation='linear')(main_input)
 output = Dense(1, activation = "linear", name = "out")(x)
 final_model = Model(inputs=[main_input], outputs=[output])
 final_model.compile(optimizer='adam', loss='mse')
return final_model

然后训练100次。

model = get_model(len(X_train[0]))
history = model.fit(X_train, Y_train,
 epochs = 100,
 batch_size = 64,
 verbose=1,
 validation_data=(X_test, Y_test),
 callbacks=[reduce_lr, checkpointer],
 shuffle=True)

得到以下结果：

30天的Keras神经网络预测

这个结果甚至比简单贝叶斯回归还糟糕，而且这个模型不能得到确定性估计，更重要的是，这个模型甚至不能正则化。

贝叶斯神经网络

现在我想在PyTorch中定义一个和我们在Keras中训练的相同的神经网络。

class Net(torch.nn.Module):
 def __init__(self, n_feature, n_hidden):
 super(Net, self).__init__()
 self.hidden = torch.nn.Linear(n_feature, n_hidden) # hidden layer
 self.predict = torch.nn.Linear(n_hidden, 1) # output layer
def forward(self, x):
 x = self.hidden(x)
 x = self.predict(x)
 return x

与贝叶斯回归模型相比，我们现在有两组参数（从输入到隐藏层以及从隐藏层到输出），所以我们稍微改变一下模型分布和先验：

priors = {'hidden.weight': w_prior,
 'hidden.bias': b_prior,
 'predict.weight': w_prior2,
 'predict.bias': b_prior2}

以及这个指导：

dists = {'hidden.weight': w_dist,
 'hidden.bias': b_dist,
 'predict.weight': w_dist2,
 'predict.bias': b_dist2}

不要忘记为模型中所有的分布设置不同的名字，因为不能有任何的歧义和重复！可以在源代码中查看更多细节。

源代码：

https://github.com/Rachnog/Deep-Trading/tree/master/bayesian

在拟合模型和采样后，让我们直接看最终结果：

30天的Pyro神经网络预测

这个结果看上去比之前的结果都要好！

考虑下从贝叶斯模型中学到的正则化或者权重的性质，与普通神经网络做比较，我还会看一下权重统计。下面是我在Pryo模型中如何检查参数的：

for name in pyro.get_param_store().get_all_param_names():
print name, pyro.param(name).data.numpy()

在Keras模型中我是这么做的：

import tensorflow as tf
sess = tf.Session()
with sess.as_default():
 tf.global_variables_initializer().run()
dense_weights, out_weights = None, None
with sess.as_default():
 for layer in model.layers:
 if len(layer.weights) > 0:
 weights = layer.get_weights()
 if 'dense' in layer.name:
 dense_weights = layer.weights[0].eval()
 if 'out' in layer.name:
 out_weights = layer.weights[0].eval()

举个例子，对于Keras模型中，最后一层的权重平均值和标准差分别是-0.0025901748，0.30395034，对于Pyro模型分别是0.0005974418和0.0005974418，数值更小，模型性能更好！

就像许多L2或者dropout这种正则化方法做的那样，让参数尽可能接近0，然后我们可以用变分推理来实现！对于隐藏层权重的情况就更有趣了。

我们把一些权重向量画出来，蓝色代表Keras的权重，橙色代表Pyro的权重：

输入和隐藏层间的一些权重

有趣的是，事实上不仅权重的均值和标准差很小，而且权重变得更加稀疏，所以基本上我们对于第一组权重用到了稀疏表示（类似L1正则），对第二组用到了类似L2正则表示，简直不可思议！不要忘记试一下代码哦！

结论

我们用了一种新奇的方法来训练神经网络。我们更新权重的分布而不是顺序更新静态权重。所以我们得到更有趣和可靠的结果。

我想要强调的是，贝叶斯方法给了我们一个机会，使得我们可以不手动添加正则项的情况下对神经网络进行正则化，理解模型的不确定性，并尽可能使用更少的数据得到更好的结果。欢迎继续关注!:)

原文链接：

https://medium.com/@alexrachnog/financial-forecasting-with-probabilistic-programming-and-pyro-db68ab1a1dba