深度学习之BP神经网络--Stata和R同步实现（附数据和代码）
凡是搞计量经济的，都关注这个号了
稿件：econometrics666@126.com
关于机器学习方法，可参阅如下文章：1机器学习方法出现在AER, JPE, QJE等顶刊上了！，2前沿: 机器学习在金融和能源经济领域的应用分类总结，3文本分析的步骤, 工具, 途径和可视化如何做？4文本大数据分析在经济学和金融学中的应用, 最全文献综述，5最全: 深度学习在经济金融管理领域的应用现状汇总与前沿瞻望, 中青年学者不能不关注！6Top前沿: 农业和应用经济学中的机器学习, 其与计量经济学的比较, 不读不懂你就out了！7机器学习和大数据计量经济学, 你必须阅读一下这篇，8机器学习与Econometrics的书籍推荐, 值得拥有的经典，9机器学习在微观计量的应用最新趋势: 大数据和因果推断，10机器学习第一书, 数据挖掘, 推理和预测，11Top, 机器学习是一种应用的计量经济学方法, 不懂将来面临淘汰危险！12最新: 运用机器学习和合成控制法研究武汉封城对空气污染和健康的影响! 13陈硕: 回顾与展望经济学研究中的机器学习，14机器学习对经济学研究的影响研究进展综述。
注：欲了解更多关于最新机器学习、深度学习等前沿计量方法和代码，可以到社群交流访问。
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，因变量）。
上图中以房价预测为例。假设我们先知道房子的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的函数图像大致如下（来源于百度百科）：
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条竖线（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函数的图像（代码省略）如下。
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）正向传播
Forward Propagation，也叫正向传播。如下面两个图所示，神经网络的正向传播就是样本信息向后传递的过程。由于Z是样本特征X的函数，A是Z的函数，相当于A是X的复合函数。同时，各隐藏层拥有本层的Z和A函数，隐藏层之间通过A（可以看成是某层的预测值）进行信息传递（当成下一层的X）。从信号（原始数据）传递的角度来看，正向传播将输入信号通过隐藏层后，作用于输出节点，经过非线性变换转化为输出信号。下图中还给出了成本函数L()的一种形式。
正向传递中的数学算法的通用表达如下：
下图展示了更具体的前向传递向量化算法实现过程。
（2）后向传播
Backward Propagation，反向传播。一次正向传播加一次反向传播，就是一次完整的iteration迭代。经过一次正向传播后，若实际输出与真实结果差距较大，则启动误差的反向传播。
如下图所示，后向传播根据前向传播计算出来的A（激活值，A=g(Z)）、真实值Y和暂存的Z（Z=Wx+b），用求导的方式求出每一层参数的梯度，以调整模型参数值，进而降低成本函数的误差（梯度下降法）。和前向传播的复合函数构造过程对应，后向传播的实现包括了复合函数从外到内的链式求（偏）导过程。
通俗来说，反向传播是将输出误差反向逐层传递，算出所有层的误差估计，即将输出误差摊分给各层所有单元，并在各层算出误差信号后对各单元的参数进行调节（误差以某种形式在各网络层进行表示）。调整正向传播中输入层和隐藏层、隐藏层和输出层对应的参数，即可推动误差随梯度方向下降。经过重复的学习，确定和最小误差所对应的参数后，学习停止。
小结：正向传播求损失，反向传播传误差，基于误差值进行参数更新；反复迭代求出最优模型。
神经网络算法的数学推导有一定的学习门槛，对于初学者来说需要一定投入才能弄懂。吴恩达在他的深度学习系列课程中提到，神经网络的反向传递的数学推断几乎是是机器学习中最复杂的，需要用到矩阵的运算、矩阵的求导和链式求导法则等数学知识。
1.4 参数和超参数
神经网络模型的参数主要包括W权重矩阵和b偏误向量，前文已提到过。除了参数，神经网络模型还包括多个超参数。超参数是决定参数大小的参数，主要包括：
Learning Rate α学习率α。如下图所示，学习率可以看成是权重w下降部分（成本函数J对w的导数）的系数，其取值介于0和1之间。一个好的学习率，可以让成本函数下降更快，收敛在更低的水平。不过也要注意，成本函数的下降速度过快也可能带来问题。若学习率取得过大，算法虽然收敛更快，但更容易陷入局部最优；反之若取值较小，算法的收敛速度较慢，但可以一步步逼近全局最优。
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 神经网络与R实现
在R语言中，实现BP神经网络建模的包和函数非常丰富，包括caret包中的nnet函数和neuralnet包中的neuralnet函数，以及本文使用的RSNNR包中的mlp函数。在mlp函数的具体实现中，我们使用了Rprop学习函数。Rprop学习函数的全称是Resilient Back Propagation, 即弹性反向传播。Rprop是普通反向传播BP算法的优化，它在两方面做出改进：
训练速度通常更快； Rprop不需要指定learning rate学习率等超参数。
3 分类：iris鸢尾花类型的预测
3.1 iris鸢尾花分类R实现
首先让我们使用之前用过的iris数据，及nnet包中的nnet函数建立一个简单的神经网络模型。模型中的输出（标签）变量是Species，其它变量为特征（自）变量。我们的目标的建立一个简单的神经网络学习算法，用来对鸢尾花进行分类。同时，我们还将该算法的性能和基于交叉验证的KNN算法进行比较。
iris <- iris
> head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
> library(nnet)
> detach("package:datasets", unload = TRUE) #这是加载nnet包后自动执行的代码
> iris.nn <- nnet(Species~., data=iris, size=2) #2个隐藏层
# weights: 19
initial value 189.255340
iter 10 value 77.905956
iter 20 value 75.777510
iter 30 value 75.294014
iter 40 value 71.919854
iter 50 value 69.534235
iter 60 value 64.106045
iter 70 value 60.466622
iter 80 value 40.041381
iter 90 value 9.411000
iter 100 value 6.262853
final value 6.262853
stopped after 100 iterations
> Species_hat <- predict(iris.nn, iris, type="class")
> table(Species_hat, iris$Species)
Species_hat setosa versicolor virginica
setosa 50 0 0
versicolor 0 49 1
virginica 0 1 49
> print( (150-2)/150) #预测准确率为98.67%。下面将使用KNN算法建模对比
[1] 0.9866667
从最后结果我们可以看出，结果交叉表中只有两个元素不为0和为2，预测准确率达到98.67%，说明该2个隐藏层的神经网络模型的预测能力非常不错。
3.2 和KNN算法对比
为了确定最优的KNN模型，考虑使用Cross Validation交叉验证搜寻最优的KNN模型（之前关于KNN算法介绍的推文中没有介绍到优化的实现，在此算是补上）。
> library(caret)
> library(ggplot2)
> set.seed(1898)
> trcl <- trainControl(method="cv", number=5) #5重交叉验证
> trgrid <- expand.grid(k=seq(1,50,2))
> irisknnfit <- train(x=iris[, -5], y=iris$Species, method="knn",
+ trControl = trcl, tuneGrid=trgrid)
> irisknnfit
k-Nearest Neighbors
150 samples
4 predictor
3 classes: 'setosa', 'versicolor', 'virginica'
No pre-processing
Resampling: Cross-Validated (5 fold)
Summary of sample sizes: 120, 120, 120, 120, 120
Resampling results across tuning parameters:
k Accuracy Kappa
1 0.9533333 0.93
3 0.9600000 0.94
5 0.9533333 0.93
7 0.9666667 0.95
9 0.9666667 0.95
11 0.9600000 0.94
13 0.9800000 0.97 #k=13时的精确度性能表现最优
15 0.9666667 0.95
17 0.9666667 0.95
19 0.9533333 0.93
21 0.9600000 0.94
23 0.9466667 0.92
25 0.9466667 0.92
27 0.9333333 0.90
29 0.9466667 0.92
31 0.9333333 0.90
33 0.9400000 0.91
35 0.9400000 0.91
37 0.9400000 0.91
39 0.9466667 0.92
41 0.9333333 0.90
43 0.9266667 0.89
45 0.9200000 0.88
47 0.9200000 0.88
49 0.9200000 0.88
Accuracy was used to select the optimal model using the largest value.
The final value used for the model was k = 13.
> plot(irisknnfit, main="KNN", family="STKaiti") #画图展示参数搜索结果
KNN算法中当k=13时模型的预测精确度最高（98.00%），但仍比不上上文中的2个隐藏层的神经网络模型的98.67%。通常而言，神经网络算法的性能是优于KNN的。
4 回归：工资BPNN模型分析和预测
先简要说明下本案例所用到的prestige数据。该数据来自R语言中的car包，包括一些和职业及职业声望有关的特征变量，以及结果变量收入。为了简化建模过程，本案例只使用其中的3个变量：首先是结果变量income，表示年收入，单位为美元；第一个自变量是education受教育年限，单位是年；第二个自变量为prestige，可当成是职业声望评分。三个变量全部为连续变量。
4.1 R神经网络建模
第一步，准备相应的包和数据。在此使用RSNNS包中的mlp函数建立BP神经网络模型。
> library(car); library(carData) #数据prestige所在的包
> library(RSNNS) #本案例的核心包，用于BP神经网络回归建模
> library(caret) #机器学期集成包
> library(NeuralNetTools) #用于神经网络画图的包
> library(Metrics) #用于计算RMSE
> prestige <- Prestige[c(1,4,2)] #只保留3列：education-prestige-income
> head(prestige) #查看数据的头部
education prestige income
gov.administrators 13.11 68.8 12351
general.managers 12.26 69.1 25879
accountants 12.77 63.4 9271
purchasing.officers 11.42 56.8 8865
chemists 14.62 73.5 8403
physicists 15.64 77.6 11030
> dim(prestige) #102行*3列
[1] 102 3
> summary(prestige$education); #查看3个变量的情况
Min. 1st Qu. Median Mean 3rd Qu. Max.
6.380 8.445 10.540 10.738 12.648 15.970
> summary(prestige$prestige);
Min. 1st Qu. Median Mean 3rd Qu. Max.
14.80 35.23 43.60 46.83 59.27 87.20
> summary(prestige$income)
Min. 1st Qu. Median Mean 3rd Qu. Max.
611 4106 5930 6798 8187 25879
> #数据的归一化处理
> #自变量education和prestige的归一化
> prestige[c(1,2)] <- normalizeData(prestige[c(1, 2)], "0_1")
> #因变量income的归一化，并提取归一化参数
> income <- normalizeData(prestige$income, "0_1")
> NormParameters <- getNormParameters(income)
> summary(prestige$education); #查看归一化后的效果
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.2153 0.4338 0.4544 0.6535 1.0000
> summary(prestige$prestige);
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.2821 0.3978 0.4424 0.6143 1.0000
> summary(income)
V1
Min. :0.0000
1st Qu.:0.1383
Median :0.2105
Mean :0.2449
3rd Qu.:0.2998
Max. :1.0000
> #使用RSNNS包中的splitForTrainingAndTest函数，将数据进行随机切割
> set.seed(1898)
> datasplist <- splitForTrainingAndTest(prestige[c(1,2)], income, ratio=0.3) #30%数据用于测试，70%数据用于训练
> dim(datasplist$inputsTrain) #训练数据的自变量，71行*2列
[1] 71 2
> dim(datasplist$inputsTest) #测试数据的31行*2列
第二步，基于测试数据，使用mlp()函数建立BP神经网络回归模型。
> BPNNreg <- mlp(datasplist$inputsTrain, datasplist$targetsTrain, maxit=200,
size=c(10, 10), learnFunc="Rprop",
inputsTest=datasplist$inputsTest,
targetsTest=datasplist$targetsTest,
metric="RMSE")
参数说明：（1）第一个和第二个参数表示训练数据的自变量（inputs）和因变量（outputs）数据；（2）maxit=200表示max iteration最大迭代次数是200（后面会画图检验）；（3）size表示共2层隐藏网络，每层有10个神经元；（4）learnFunc表示神经网络的学习函数类型；（5）第六和第七个参数表示测试数据的自变量和因变量数据；（6）最后一个参数metric=”RMSE”表示衡量误差的方法是均方误差。其它省略的参数使用的是默认值。
第三步，查看已训练好的神经网络回归模型。
Class: mlp->rsnns
Number of inputs: 2
Number of outputs: 1
Maximal iterations: 200
Initialization function: Randomize_Weights
Initialization function parameters: -0.3 0.3
Learning function: Rprop
Learning function parameters: 0.2 0
Update function:Topological_Order
Update function parameters: 0
Patterns are shuffled internally: TRUE
Compute error in every iteration: TRUE
Architecture Parameters:
$size
[1] 10 10
All members of model:
[1] "nInputs" "maxit" "initFunc" "initFuncParams" "learnFunc" "learnFuncParams"
[7] "updateFunc" "updateFuncParams" "shufflePatterns" "computeIterativeError" "snnsObject" "archParams"
[13] "IterativeFitError" "IterativeTestError" "fitted.values" "fittedTestValues" "nOutputs"
> summary(BPNNreg) #summary函数会展示更多的模型细节。在此省略其结果
> par(cex=0.6)
> plotnet(BPNNreg, pos_col="red", neg_col="grey")
# 这里pos_col表示正的连接权重的颜色是红色，负的为灰色
# 此外，在层数和神经节点较多的情况下，画图很可能比建模更费时 > #用plotIterativeError函数画出模型的迭代优化过程
> plotIterativeError(BPNNreg, main="BPNN回归迭代误差") 图中黑线表示训练集的迭代误差，迭代次数在150次左右达到稳定；红线表示测试集的迭代误差，在10次以内迭代误差就达到了稳定。这说明之前建模时候设置maxit=200是合理的。
> #使用plotRegressionError()函数，可视化模型的拟合误差
> #训练集的拟合误差
> plotRegressionError(datasplist$targetsTrain, BPNNreg$fitted.values, main="BPNN regssion train fit")
图中红色拟合线越接近黑色的对角线，表示拟合效果越好。
> #测试集的拟合误差
> plotRegressionError(datasplist$targetsTest, BPNNreg$fittedTestValues, main="BPNN regssion test fit")
第四步，使用测试数据检验BP神经网络模型的均方误差大小。
train_rmse <- Metrics::rmse(BPNNregpre_train, BPNN_train) #训练数据集上的均方误差
> train_rmse #3131.223
[1] 3131.223
> #查看测试集上的RMSE
> BPNNregpre_test <- denormalizeData(BPNNreg$fittedTestValues, NormParameters) #测试数据中拟合的y还原
> BPNN_test <- denormalizeData(datasplist$targetsTest, NormParameters) #测试数据中的y还原
> test_rmse <- Metrics::rmse(BPNNregpre_test, BPNN_test) #训练数据的均方误
> test_rmse #1736.741
[1] 1736.741
4.2 Stata神经网络建模
第一步，数据准备。同样使用的是prestige数据，但随机切割后一般会有些许差别。
. clear
. set seed 1898
. use "D:\R\prestige.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:\R\BPNNm1.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
第四步，查看神经网络模型中自变量的伪边际效应
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)
*从画图结果可以看出，随着受教育程度的增加，收入整体是上升的
第五步，准备测试数据。
*数据准备
. clear
. set seed 1898
. use "D:\R\prestige.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:\R\BPNNm1.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:\R\prestige.train.dta", clear
. reg income prestige education //RMSE=3054.5
. use "D:\R\prestige.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
4.4 结果汇总
| 工具-模型 | 训练模型的RMSE | 预测的RMSE |
|---|---|---|
| Stata-OLS | 3054.5 | 2856.00 |
| Stata-BP神经网络 | NA | 3400.33 |
| R-Rprop神经网络 | 3131.22 | 1736.74 |
从上表可以大致看出，R中的Rprop弹性BP神经网络的预测效果最好，和预期一致。不过也需注意到Rprop神经网络训练模型的RMSE比OLS的更大。理论上来说，不管是训练数据还是预测数据，Rprop模型的性能应该最优，Stata中的BP神经网络次之，OLS最后。
5 小结
和其它一些机器学习模型一样，寻找好的神经网络模型的过程依靠的往往是经验法则，需要经过反复的参数和超参数的调整和比较。最优超参数的选择，还可能随着应用领域、数据更新和硬件技术的升级而改变。当然，短期内的模型优化方法主要还是Cross Validation交叉验证，用现有数据重复地检验模型的性能表现。本案例中展示的是一个较为简单的BP神经网络，主要目的是介绍相关的建模过程。在掌握了基本的建模过程后，大家可以尝试从以下几个方面建立更复杂的模型：
引入更多自变量和观测值。本文附带有一个名为“ house_1.05m.csv”的真实房价数据拥有13个变量和约105万个观测值，可供大家进一步练习。引入 因子型（分类）变量。在引入连续变量的时候需要去量纲化（0-1或者标准化），而引入因子型变量的时候处理过程稍微更复杂，需要将分类变量转换为一系列的哑变量后纳入回归模型。该转换可以通过RSNNS包的decodeClassLabels函数实现。将隐藏层参数设置为 更多的层数及更多的神经元数（可能需要换一台性能更好的电脑）。在此顺便向大家展示下prestige回归案例中c(50, 50, 50, 50)的神经网络长什么样。
参考文献
陈强，《高级计量经济学》，第二版 薛震，孙玉林，《R语言统计分析与机器学习》，第一版 Scott V. Burger，《基于R语言的机器学习》，第一版 Robert I. Kabacoff，《R in Action》，第一版 Rumelhart, D., Hinton, G. & Williams, R. Learning representations by back-propagating errors. Nature 323, 533–536 (1986). https://doi.org/10.1038/323533a0 Thorsten Doherr, 2018. "BRAIN: Stata module to provide neural network," Statistical Software Components S458566, Boston College Department of Economics. Thorsten Doherr, codes of Stata command “brain”, https://github.com/ThorstenDoherr/brain James McCaffrey, 2015, How To Use Resilient Back Propagation To Train Neural Networks, https://visualstudiomagazine.com/articles/2015/03/01/resilient-back-propagation.aspx Riedmiller M. (1994) Rprop - Description and Implementation Details. Technical Report. University of Karlsruhe. Stata中“brain”命令的帮助文件。 吴恩达的深度学习系列课程，在Coursera、网易云课堂和B站均可公开访问。 R中mlp()、nnet()和neuralnet()等函数的帮助文件。
以上为数据和Stata及R代码。
关于相关计量方法视频课程，文章，数据和代码，参看 1.面板数据方法免费课程, 文章, 数据和代码全在这里, 优秀学人好好收藏学习！2.双重差分DID方法免费课程, 文章, 数据和代码全在这里, 优秀学人必须收藏学习！3.工具变量IV估计免费课程, 文章, 数据和代码全在这里, 不学习可不要后悔！4.各种匹配方法免费课程, 文章, 数据和代码全在这里, 掌握匹配方法不是梦！5.断点回归RD和合成控制法SCM免费课程, 文章, 数据和代码全在这里, 有必要认真研究学习！6.空间计量免费课程, 文章, 数据和代码全在这里, 空间相关学者注意查收！7.Stata, R和Python视频课程, 文章, 数据和代码全在这里, 真的受用无穷！
2.5年，计量经济圈近1000篇不重类计量文章，
可直接在公众号菜单栏搜索任何计量相关问题,
Econometrics Circle
数据系列：空间矩阵 | 工企数据 | PM2.5 | 市场化指数 | CO2数据 | 夜间灯光 | 官员方言 | 微观数据 | 内部数据计量系列：匹配方法 | 内生性 | 工具变量 | DID | 面板数据 | 常用TOOL | 中介调节 | 时间序列 | RDD断点 | 合成控制 | 200篇合辑 | 因果识别 | 社会网络 | 空间DID数据处理：Stata | R | Python | 缺失值 | CHIP/ CHNS/CHARLS/CFPS/CGSS等 |干货系列：能源环境 | 效率研究 | 空间计量 | 国际经贸 | 计量软件 | 商科研究 | 机器学习 | SSCI | CSSCI | SSCI查询 | 名家经验计量经济圈组织了一个计量社群，有如下特征：热情互助最多、前沿趋势最多、社科资料最多、社科数据最多、科研牛人最多、海外名校最多。因此，建议积极进取和有强烈研习激情的中青年学者到社群交流探讨，始终坚信优秀是通过感染优秀而互相成就彼此的。