盲区行者：深度学习之BP神经网络--Stata和R同步实现（附R数据和代码）zhuanlan.zhihu.com

原公众号推文标题：深度学习之BP神经网络-Stata和R同步实现（附数据和代码）

神经网络（Neural Network，或Artificial Neural Network，简称NN或ANN）是Deep Learning深度学习（DL属于ML，ML属于AI）算法的代表。和SVM支持向量机类似，神经网络属于典型的black box黑箱方法。从更广泛的角度来看，神经网络可以实现几乎所有的学习算法——既可以应用于有监督学习的分类和回归，也可以用于无监督学习的聚类等。

神经网络的分类有多种，从演进的角度来看，大概包括早期的单层神经网络（也叫感知器，Perceptron）和前馈神经网络（Feedforward Neural Network，FNN），到BP（Backward Propagation，反向传播）神经网络，再到现在更为复杂的卷帙（Convolutional）、循环（Recurrent）、深度信念（Deep Belief）和生成对抗（Generative Adversarial）神经网络等等。

早期的前馈神经网络由于过于简单，应用价值不大。随后的BP算法的提出，使得神经网络的性能得到了大幅提升。在此基础之上演化出来了多种神经网络算法，大规模应用于Computer Vision计算机视觉、Speech Recognition语音识别和Natural Language Processing自然语言处理，以及在线广告、网页搜索和商品推荐等。

BP神经网络是神经网络的典型代表，是初学者的重点学习对象。此外，Stata中目前能实现神经网络建模的的只有BP神经网络brain命令，尚无其它神经网络算法的引进，故本文在此重点介绍BP神经网络。

1 神经网络介绍

1.1 神经网络原理

BP神经网络属于监督学习，由激活函数、网络结构和优化算法三大部分组成。对于初学者来说，我们也可以把BP神经网络模型大概当成一个庞杂的复合函数。

BP神经网络是在前馈FNN神经网络模型的基础上改进。下面的视频截图（本文中引用的视频截图全部来自吴恩达深度学习课程）展示的是单隐藏层的神经网络，便于我们大概了解FNN的原理。这里所谓前馈，是指向前传播。信息由输入变量（inputs，自变量）开始，往前经过中间的各个隐藏层向前传导，最终到达输出节点（outputs，因变量）。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_R语言 BP神经网络时间序列预测

上图中以房价预测为例。假设我们先知道房子的size面积、bedroom卧室的数量、zipcode邮编和当地的wealth人均收入，打算基于这4个自变量来预测房价price。这里自变量中的size和bedroom的数量，可能会通过购买者的family size来影响房价；邮编zipcode可能通过该社区的walkability便利性来影响房价；zipcode还可能和wealth一起，通过当地的school quality影响房价。这个推导过程和计量经济学中的中介变量、结构化方程SEM和RCT中的因果链推导原理有类似之处。

在完成前馈（正向）传播之后，信息再由输出节点反向传播。在介绍正向和反向传播具体原理之前，我们先简单了解下常见的三种激活函数。

1.2 三种常见的激活函数

Activation Function激活函数用于“在神经网络模型中加入非线性因素，以解决线性模型所不能解决的问题”。激活函数也用于反向传递中的参数修正，所以通常要求是可微的。激活函数的形式可以由用户自行定义，而常见的主要是以下三个。

（1）ReLU修正线性单元

ReLU，全称是Rectified Linear Unit，是最基础的一种激活函数。相对于后面两种激活函数，ReLU最大的优点是斜率的稳定性。不管x的取值多大，在梯度下降法中的下降斜率都保持一个较大的值，而不像Sigmoid或Tanh函数一样可能趋近于零而导致梯度算法下降较小，最终算法学习的速度会变得很慢。 ReLU的缺点也较明显，即当自变量取值小于0的时候其函数斜率恒等于0。因此，ReLU在自变量小于0的区域的斜率被改进成为了某个较小的常数，比如0.01。这种改进后的激活函数被称为Leaky ReLU，泄露ReLU。ReLU的函数图像大致如下（来源于百度百科）：

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_02

ReLU和Leaky ReLU通常用于卷积神经网络和循环神经网络等。

（2）Sigmoid函数

机器学习的Sigmoid函数，又被称为Logit函数，或Logistics函数，在计量中常用于结果变量是0-1等分类变量的建模。其函数形式如下：

y= 1/(1- e^(-x) )

接下来我们用Stata画出该函数的图形，让大家有个直观的了解。

. clear
. set obs 1000000 
. gen x = runiform(-20, 20)
. gen y = 1/(1+exp(-x))
. line y x, sort xline(-5) xline(0) xline(5) xlabel(-20(5)20)
*在x=-5/0/5初加了3条竖线

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络预测python代码_03

3）Tanh双曲正切函数

Tanh函数形式如下：

y = sinhx/coshx=((e^x-e^(-x))/((e^x+e^(-x))

Tanh函数可以看成是Sigmoid函数的改进，通常只要能用Sigmoid函数的地方，Tanh函数的表现都会更优。当然，由于Sigmoid的0-1的结果可能，让它在0-1分类的算法中无可替代。在具体操作中，涉及0-1分类的神经网络通常在输出层使用Sigmoid函数，其它层使用Tanh或其它激活函数。 在此也用Stata画出Tanh函数的图像（代码省略）如下。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_04

1.3正向传播和反向传播

先简要说明涉及到的函数和数学符号：数学标记符号中大写表示向量，如x表示某个具体的自变量，而X表示自变量矩阵；Z=WX+b，这里W指weight权重矩阵，b表示bias偏误向量；σ()或g()表示激活函数，σ(Z)=A（A表示activation，输出层的A=Yhat），信息正是通过Z=WX+b和σ(Z)=A逐层向后传播；L(A, Y)中的L表示loss，为成本（基于损失函数）函数，确定Y和Yhat之间的误差。

正向传播和反向传播再具体介绍BP神经网络的模型算法。

（1）正向传播

正向传播。如下面两个图所示，神经网络的正向传播就是样本信息向后传递的过程。由于Z是样本特征X的函数，A是Z的函数，相当于A是X的复合函数。同时，各隐藏层拥有本层的Z和A函数，隐藏层之间通过A（可以看成是某层的预测值）进行信息传递（当成下一层的X）。从信号（原始数据）传递的角度来看，正向传播将输入信号通过隐藏层后，作用于输出节点，经过非线性变换转化

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_05

正向传递中的数学算法的通用表达如下：

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_06

下图展示了更具体的前向传递向量化算法实现过程。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络预测python代码_07

（2）后向传播

Backward Propagation，反向传播。一次正向传播加一次反向传播，就是一次完整的iteration迭代。经过一次正向传播后，若实际输出与真实结果差距较大，则启动误差的反向传播。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_R语言 BP神经网络时间序列预测_08

如下图所示，后向传播根据前向传播计算出来的A（激活值，A=g(Z)）、真实值Y和暂存的Z（Z=Wx+b），用求导的方式求出每一层参数的梯度，以调整模型参数值，进而降低成本函数的误差（梯度下降法）。和前向传播的复合函数构造过程对应，后向传播的实现包括了复合函数从外到内的链式求（偏）导过程。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_09

通俗来说，反向传播是将输出误差反向逐层传递，算出所有层的误差估计，即将输出误差摊分给各层所有单元，并在各层算出误差信号后对各单元的参数进行调节（误差以某种形式在各网络层进行表示）。调整正向传播中输入层和隐藏层、隐藏层和输出层对应的参数，即可推动误差随梯度方向下降。经过重复的学习，确定和最小误差所对应的参数后，学习停止。 　

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_10

小结：正向传播求损失，反向传播传误差，基于误差值进行参数更新；反复迭代求出最优模型。

神经网络算法的数学推导有一定的学习门槛，对于初学者来说需要一定投入才能弄懂。吴恩达在他的深度学习系列课程中提到，神经网络的反向传递的数学推断几乎是是机器学习中最复杂的，需要用到矩阵的运算、矩阵的求导和链式求导法则等数学知识。

1.4 参数和超参数

神经网络模型的参数主要包括W权重矩阵和b偏误向量，前文已提到过。除了参数，神经网络模型还包括多个超参数。超参数是决定参数大小的参数，主要包括：

  • Learning Rate α 学习率α。如下图所示，学习率可以看成是权重w下降部分（成本函数J对w的导数）的系数，其取值介于0和1之间。一个好的学习率，可以让成本函数下降更快，收敛在更低的水平。不过也要注意，成本函数的下降速度过快也可能带来问题。若学习率取得过大，算法虽然收敛更快，但更容易陷入局部最优；反之若取值较小，算法的收敛速度较慢，但可以一步步逼近全局最优。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_R语言 BP神经网络时间序列预测_11

  • Training Factor训练因子。在Stata的brain命令中叫“eta”（默认值是0.25），用来指定学习率的上限和下限。
  • 隐藏层层数，the number of hidden layer。
  • Size，某个隐藏层神经元的个数。

2 神经网络的实现

2.1 神经网络与Stata实现

目前Stata 16中用于神经网络建模的是非官方命令“brain”，由Thorsten Doherr于2018年通过SSC和Github发布。该命令基于Stata的Mata矩阵编程语言，建立backward propagation向后传播型的multi-layer多层神经网络。Backwark propagation简称BP，故该算法也被称为BP神经网络。BP神经网络最早在1986年由Rumelhart和McCelland发表在Nature上的论文《Learning Representations by Back-propagating Errors》提出。

Doherr对brain命令的简介：“ brain”为实现反向传播神经网络算法而设计，简洁易上手。训练后的神经网络模型可以通过.brn文件保存或调用。具体模型结果以公开且可访问的矩阵的形式展示，方便老版本Stata的读取。 brain命令还附带很多函数，以便于 pseudo-marginal effects 伪边际效应的分析；但它的主要功能还是在于预测，即倾向得分计算或者分类。和其它机器学习算法类似，边际效应的计算，是因果推断中讨论自变量对因变量影响程度的关键。

2.2 Stata神经网络建模

第一步，数据准备。同样使用的是prestige数据，但随机切割后一般会有些许差别。

clear
. set seed 1898
. use "D:Rprestige.train.dta", clear

. sum prestige 
. dis r(min) //14.8	
. dis (r(max)-r(min)) //72.399997
. replace prestige=(prestige-r(min))/(r(max)-r(min)) //prestige归一化
. sum prestige

. sum education 
. dis r(min) //6.3800001
. dis r(max)-r(min) //9.5900002
. replace education=(education-r(min))/(r(max)-r(min)) //education归一化
. sum education

. sum income 
. dis r(min) //918
. dis r(max)-r(min) //24961
. replace income = (income-r(min)) / (r(max) - r(min)) //income归一化
. sum income

第二步，建立BP神经网络模型。

. clear
*定义BP神经网络模型
. brain define, input(prestige education) output(income) hidden(10 10)
Defined matrices:
   input[4,2]
  output[4,1]
  neuron[1,23]
   layer[1,4]
   brain[1,151] 

*输入变量是自变量prestige和education，输出变量是income
*hidden(10 10)表示共2个隐藏层，每个各10个神经节点。和R部分一致

. brain train, iter(200) eta(0.25) nosort 
*iter表示迭代的次数，在此设置200次
*eta为训练因子，默认取值为0.25，表示学习率的初始值均匀分布于[-0.25, 0.25]
*nosort表示在训练前不对数据进行排序

. brain save "D:RBPNNm1.brn"  //保存训练好的BP神经网络，方便后面调用进行预测

第三步，查看已经建立好的BP神经网络模型。

. matrix layer //查看各层神经节点数量
 layer[1,4]
            in  hid1  hid2   out
 neurons     2    10    10     1 */

. matrix list input  //查看输入变量
input[4,2]
         prestige  education
   min          0          0
  norm          1          1
 value          0          0
signal          0          0 */

. matrix list output //查看输出变量
       income
   min       0
  norm       1
 value       0
signal       0 

. matrix list neuron //查看神经节点
neuron[1,23]
          in1    in2   h1n1   h1n2   h1n3   h1n4   h1n5   h1n6   h1n7   h1n8   h1n9  h1n10   h2n1
signal      0      0      0      0      0      0      0      0      0      0      0      0      0

         h2n2   h2n3   h2n4   h2n5   h2n6   h2n7   h2n8   h2n9  h2n10   out1
signal      0      0      0      0      0      0      0      0      0      0 

. matrix list brain //查看各节点的连接权重
*只展示部分结果
brain[1,151]
            h1n1w1      h1n1w2       h1n1b      h1n2w1      h1n2w2       h1n2b      h1n3w1
weight   .41241582   -.4473998    .1887836   1.0385614    .4644028  -.32226553   .71407534

            h1n3w2       h1n3b      h1n4w1      h1n4w2       h1n4b      h1n5w1      h1n5w2
weight   .85327234  -.18876532  -.52128979  -.39156865   -.2958014   .88854821   .79453232

第四步，查看神经网络模型中自变量的伪边际效应

.brain margin //查看边际效应
                income
 prestige  0.074326136  
education  0.056956779  
/*结果的解读：（1）职业声望每增加一单位（从0变到1），收入增加0.0743个单位(约合1855.25美元)；（2）受教育程度每增加一单位，收入增加0.0570个单位（约合1421.6982美元）。*/
                                                                        

/*                   income
   prestige  14600.733767469
  education  -5672.844525771 */ // 未归一化的结果。education系数居然是负的！这显然和常识不符。让我们通过画图检查下数据

. twoway (line income education, sort)  (lfit income education)
*从画图结果可以看出，随着受教育程度的增加，收入整体是上升的

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络matlab代码_12

第五步，准备测试数据。

*数据准备
. clear
. set seed 1898
. use "D:Rprestige.test.dta", clear
. sum prestige 
. dis r(min) //23.200001	
. dis (r(max)-r(min)) //59.100002
. replace prestige = (prestige-r(min)) / (r(max) - r(min)) //prestige归一化处理
. sum prestige

. sum education 
/*
    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
   education |         28    10.67536    2.965029        6.6      15.94 */
   
. dis r(min) //6.5999999
. dis r(max)-r(min) //9.3399997
. replace education = (education-r(min)) / (r(max) - r(min)) //education归一化处理
. sum education

. sum income 
. dis r(min) //611
. dis r(max)-r(min) //18652
. replace income = (income-r(min)) / (r(max) - r(min)) //income归一化处理
. sum income

*导入之前训练好的神经网络
. brain load "D:RBPNNm1.brn"
Loaded matrices:
   input[4,2]
  output[4,1]
  neuron[1,23]
   layer[1,4]
   brain[1,151] 

. brain think income_hat //基于预测数据进行预测，生成income_hat变量

*还原income和income_hat
. sum income income_hat
. replace income = income*18652 + 611
. replace income_hat = income_hat*18652 + 611
. sum income income_hat
/*
    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
      income |         28    6668.821     3984.65        611      19263
  income_hat |         28    5018.372    1651.003   3066.421   8184.359 */

第六步，评价Stata BP神经网络模型预测的表现。

. gen income_se = (income_hat - income)^2
. sum income_se
. scalar income_se_mean = r(mean)
. dis sqrt(income_se_mean)  //RMSE = 3400.3261

4.3 Stata OLS建模对比

. set seed 1898
. use "D:Rprestige.train.dta", clear
. reg income prestige education //RMSE=3054.5

. use "D:Rprestige.test.dta", clear
. predict income_hat
. gen income_se = (income_hat - income)^2
. sum income_se
. scalar income_se_mean = r(mean)
. dis sqrt(income_se_mean)  //RMSE = 2855.9972

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_bp神经网络代码_13

3 小结

和其它一些机器学习模型一样，寻找好的神经网络模型的过程依靠的往往是经验法则，需要经过反复的参数和超参数的调整和比较。最优超参数的选择，还可能随着应用领域、数据更新和硬件技术的升级而改变。当然，短期内的模型优化方法主要还是Cross Validation交叉验证，用现有数据重复地检验模型的性能表现。 本案例中展示的是一个较为简单的BP神经网络，主要目的是介绍相关的建模过程。在掌握了基本的建模过程后，大家可以尝试从以下几个方面建立更复杂的模型：

  • 引入更多自变量和观测值。本文附带有一个名为“house_1.05m.csv”的真实房价数据拥有13个变量和约105万个观测值，可供大家进一步练习。
  • 引入因子型（分类）变量。在引入连续变量的时候需要去量纲化（0-1或者标准化），而引入因子型变量的时候处理过程稍微更复杂，需要将分类变量转换为一系列的哑变量后纳入回归模型。该转换可以通过RSNNS包的decodeClassLabels函数实现。
  • 将隐藏层参数设置为更多的层数及更多的神经元数（可能需要换一台性能更好的电脑）。在此顺便向大家展示下prestige回归案例中c(50, 50, 50, 50)的神经网络长什么样。

R语言 BP神经网络时间序列预测 r语言bp神经网络实例_R语言 BP神经网络时间序列预测_14

prestige回归案例中c(50, 50, 50, 50)的神经网络视觉化（基于R语言）

Ref

1. 陈强，《高级计量经济学》，第二版。

2. 薛震，孙玉林，《R语言统计分析与机器学习》，第一版。

3. Scott V. Burger，《基于R语言的机器学习》，第一版。

4. Robert I. Kabacoff，《R in Action》，第一版。

5. Rumelhart, D., Hinton, G. & Williams, R. Learning representations by back-propagating errors. Nature 323, 533–536 (1986). https://doi.org/10.1038/323533a0.

6. Thorsten Doherr, 2018. "BRAIN: Stata module to provide neural network," Statistical Software Components S458566, Boston College Department of Economics.

7. Thorsten Doherr, codes of Stata command “brain”, https://github.com/ThorstenDoherr/brain.

8. James McCaffrey, 2015, How To Use Resilient Back Propagation To Train Neural Networks, https://visualstudiomagazine.com/articles/2015/03/01/resilient-back-propagation.asp.

9. Riedmiller M. (1994) Rprop - Description and Implementation Details.

10. Stata中“brain”命令的帮助文件。

11. 吴恩达的深度学习系列课程，在Coursera、网易云课堂和B站均可公开访问。

12. R中mlp()、nnet()和neuralnet()等函数的帮助文件。

-------全文结束-------------------------------------------------------------------------------