算法五：神经网路（优化算法）
人工神经网络（ANN），简称神经网路，是一种模仿生物神经网路的结构和功能的物理模型或估算模型。神经网路由大量的人工神经元连结进行估算。大多数情况下人工神经网络能在外界信息的基础上改变内部结构，是一种自适应系统。现代神经网路是一种非线性统计性数据建模工具，常用来对输入和输出间复杂的关系进行建模，或拿来探求数据的模式。
人工神经网络从以下四个方面去模拟人的智能行为:
算法训练步骤：
1、定义变量与参数x（输入向量）,w（权值向量）,b（偏置）,y（实际输出）,d（期望输出）,a（学习率参数）
2、初始化，n=0,w=0
3、输入训练样本，对每位训练样本指定其期望输出：A类记为1，B类记为-1
4、计算实际输出y=sign(w*x+b)
5、更新残差向量w(n+1)=w(n)+a[d-y(n)]*x(n),0
6、判断，若满足收敛条件，算法结束，否则返回3
注意，其中学习率a为了残差的稳定性不应过大，为了彰显偏差对残差的修正不应过小，说到底，这是个经验问题。
从后面的表述来看，感知器对于线性可分的事例是一定收敛的，对于不可分问题，它无法实现正确分类。这里与我们上面提到的支持向量机的看法非常的相仿，只是确定分类直线的办法有所不同。可以这么说，对于线性可分的事例，支持向量机找到了“最优的”那条分类直线，而双层感知器找到了一条可行的直线。
我们以鸢尾花数据集为例，由于双层感知器是一个二分类器，所以我们将鸢尾花数据也分为两类，“setosa”与“versicolor”（将后两类均看做第2类），那么数据根据特点：花瓣宽度与长度做分类。
运行下边的代码：
#感知器训练结果：
a<-0.2
w<-rep(0,3)
iris1<-t(as.matrix(iris[,3:4]))
d<-c(rep(0,50),rep(1,100))
e<-rep(0,150)
p<-rbind(rep(1,150),iris1)
max<-100000
eps<-rep(0,100000)
i<-0
repeat{
v<-w%*%p;
y<-ifelse(sign(v)>=0,1,0);
e<-d-y;
eps[i+1]<-sum(abs(e))/length(e)
if(eps[i+1]<0.01){
print("finish:");
print(w);
break;
}
w<-w+a*(d-y)%*%t(p);
i<-i+1;
if(i>max){
print("max time loop");
print(eps[i])
print(y);
break;
}
}
#绘图程序
plot(Petal.Length~Petal.Width,xlim=c(0,3),ylim=c(0,8),
data=iris[iris$Species=="virginica",])
data1<-iris[iris$Species=="versicolor",]
points(data1$Petal.Width,data1$Petal.Length,col=2)
data2<-iris[iris$Species=="setosa",]
points(data2$Petal.Width,data2$Petal.Length,col=3)
x<-seq(0,3,0.01)
y<-x*(-w[3]/w[2])-w[1]/w[2]
lines(x,y,col=4)
#绘制每次迭代的平均绝对误差
plot(1:i,eps[1:i],type="o")
分类结果如图：
这是运行了7次得到的结果。与我们上面的支持向量机相比，显然神经网路的双层感知器分类不是这么的可信，有些弱。（原来的代码和图有一点小问题，感谢@xuchaohit的见谅）
我们可以尝试来做交叉验证，可以发觉交叉验证结果并不理想。
二、线性神经网路
尽管当训练样例线性可分时，感知器法则可以成功地找到一个权向量，但若果样例不是线性可分时它将不能收敛。因此，人们设计了另一个训练法则来克服这个不足，称为delta法则。
如果训练样本不是线性可分的，那么delta法则会收敛到目标概念的最佳近似。
delta法则的关键思想是使用梯度增长来搜索可能权向量的假定空间，以找到最佳拟合训练样例的权向量。
我们将算法描述如下：
1、定义变量与参数。x（输入向量）,w（权值向量）,b（偏置）,y（实际输出）,d（期望输出）,a（学习率参数）（为表述简便，我们可以将偏置划归残差向量中）
2、初始化w=0
3、输入样本，计算实际输出与偏差。e(n)=d-x*w(n)
4、调整残差向量w(n+1)=w(n)+a*x*e(n)
5、判断是否收敛，收敛结束，否则返回3
Hayjin证明，只要学习率a
我们还是以里面的鸢尾花数据为例来说这个问题。运行代码：
p<-rbind(rep(1,150),iris1)
d<-c(rep(0,50),rep(1,100))
w<-rep(0,3)
a<-1/max(eigen(t(p)%*%p)$values)
max<-1000
e<-rep(0,150)
eps<-rep(0,1000)
i<-0
for(i in 1:max){
v<-w%*%p;
y<-v;
e<-d-y;
eps[i+1]<-sum(e^2)/length(e)
w<-w+a*(d-y)%*%t(p);
if(i==max)
print(w)
}
得到分类直线：
相比感知器分类而言早已好了太多了，究其原因不外乎传递函数由二值阀值函数变为了线性函数，这也就是我们上面提及的delta法则会收敛到目标概念的最佳近似。增量法则渐近收敛到最小偏差假定，可能须要无限的时间，但无论训练样例是否线性可分就会收敛。
为了明了这一点我们考虑鸢尾花数据后两类花的分类（这里我们将前两类看做一类），使用感知器：
使用线性分类器：
但是要解释的一点是，收敛并不意味着分类疗效更好，要解决线性不可分问题须要的是添加非线性输入或则降低神经元。我们以Minsky & Papert (1969)提出的异或事例为例说明这一点。
使用线性神经网路，代码与前面完全相同，略。
第一个神经元输出：
权值： [,1] [,2] [,3]
[1,] 0.750.5 -0.5
测试： [,1] [,2] [,3][,4]
[1,] 1 0 1 1
第二个神经元输出：
权值： [,1] [,2] [,3]
[1,] 0.75-0.5 0.5
测试： [,1] [,2] [,3][,4]
[1,] 1 1 0 1
求解异或逻辑（相同取0，不同取1）有结果：（代码xor(c(1,0,1,1),c(1,1,0,1))）
[1] FALSE TRUETRUE FALSE
即0，1，1，0，分类正确。
最后再说一点，Delta规则只能训练双层网路，但这不会对其功能引起很大的影响。从理论上说，多层神经网路并不比双层神经网路更强大，他们具有同样的能力。
三、BP神经网络1、sigmoid函数分类
回顾我们上面提及的感知器，它使用示性函数作为分类的办法。然而示性函数作为分类器它的跳点让人感觉很难处理，幸好sigmoid函数y=1/(1+e^-x)有类似的性质，且有着光滑性这一优良性质。我们通过右图可以看到sigmoid函数的图象：
Sigmoid函数有着估算代价不高，易于理解与实现的优点但也有着欠拟合，分类精度不高的特点人工神经网络和支持向量机优点，我们在支持向量机一章中就可以看见sigmoid函数差劲的分类结果。
2、BP神经网络结构
BP (Back Propagation)神经网路，即偏差反传误差反向传播算法的学习过程，由信息的正向传播和偏差的反向传播两个过程组成。由右图可知，BP神经网络是一个三层的网路:
当实际输出与期望输出不符时，进入偏差的反向传播阶段。误差通过输出层，按偏差梯度增长的形式修正各层残差，向隐藏层、输入层逐层反传。周而复始的信息正向传播和偏差反向传播过程，是各层残差不断调整的过程，也是神经网路学习训练的过程，此过程仍然进行到网路输出的偏差降低到可以接受的程度，或者预先设定的学习次数为止。
3、反向传播算法
反向传播这一算法把我们上面提及的delta规则的剖析扩充到了带有隐藏节点的神经网路。为了理解这个问题，设想Bob给Alice讲了一个故事，然后Alice又讲给了Ted，Ted检测了这个事实真相，发现这个故事是错误的。现在 Ted 需要找出什么错误是Bob导致的而什么又迁怒于Alice。当输出节点从隐藏节点获得输入，网络发觉出现了偏差，权系数的调整须要一个算法来找出整个偏差是由多少不同的节点导致的，网络须要问，“是谁让我误入歧途？到如何的程度？如何填补？”这时，网络该如何做呢？
同样始于梯度降落原理，在权系数调整剖析中的惟一不同是涉及到t(p,n)与y(p,n)的差分。通常来说Wi的改变在于：
alpha * s'(a(p,n)) * d(n) *X(p,i,n)
其中d(n)是隐藏节点n的函数，让我们来看：
一方面，n 影响一个输出节点越多，n 造成网路整体的偏差也越多。另一方面，如果输出节点影响网路整体的偏差越少，n 对输出节点的影响也相应减小。这里d(j)是对网路的整体偏差的基值，W(n,j) 是 n 对 j 造成的影响，d(j) * W(n,j) 是这两种影响的总和。但是 n 几乎总是影响多个输出节点，也许会影响每一个输出结点，这样，d(n) 可以表示为：SUM(d(j)*W(n,j))
这里j是一个从n获得输入的输出节点，联系上去，我们就得到了一个培训规则。
这里每位从n接收输入的输出节点j都不同。关于反向传播算法的基本情况大致这么。
通常把第 1部份称为正向传播，把第2部份称为反向传播。反向传播的名子由此而至。
4、最速增长法与其改进
最速下降法的基本思想是：要找到某函数的最小值，最好的办法是沿函数的梯度方向探访，如果梯度记为d,那么迭代公式可写为w=w-alpha*d，其中alpha可理解为我们上面提及的学习速度。
最速下降法有着收敛速率慢（因为每次搜索与前一次均正交，收敛是锯齿形的），容易深陷局部最小值等缺点，所以他的改进办法也有不少，最常见的是降低动量项与学习率可变。
增加冲量项（Momentum）
改变学习率
5、BP神经网络的实现
（1）数据读入，这里我们还是使用R的外置数据——鸢尾花数据，由于神经网路本质是2分类的，所以我们将鸢尾花数据也分为两类（将前两类均看做第2类），按照特点：花瓣宽度与长度做分类。
（2）划分训练数据与测试数据
（3）初始化BP网路，采用包含一个蕴涵层的神经网路，训练方式使用包含动量的最速下降法，传递函数使用sigmoid函数。
（4）输入样本，对样本进行归一化，计算偏差，求解误差平方和
（5）判断是否收敛
（6）根据偏差调整残差。权值按照以下公式进行调整：
Delta（w）= alpha *s'(a(p,n))*(t(p,n) - y(p,n)) * X(p,n,j)
其中，alpha为学习率，s'(a(p,n))*(t(p,n)- y(p,n))为局部梯度。此外，由于使用了有动量因子的最速下降法，除第一次外，后续改变量应为：
Delta（w）(n)=-alpha*(1-mc)*Delta(w)(n)+mc*Delta(w)(n-1)
(7)测试，输出分类正确率。
完整的R代码：
iris1<-as.matrix(iris[,3:4])
iris1<-cbind(iris1,c(rep(1,100),rep(0,50)))
set.seed(5)
n<-length(iris1[,1])
samp<-sample(1:n,n/5)
traind<-iris1[-samp,c(1,2)]
train1<-iris1[-samp,3]
testd<-iris1[samp,c(1,2)]
test1<-iris1[samp,3]
set.seed(1)
ntrainnum<-120
nsampdim<-2
net.nin<-2
net.nhidden<-3
net.nout<-1
w<-2*matrix(runif(net.nhidden*net.nin)-0.5,net.nhidden,net.nin)
b<-2*(runif(net.nhidden)-0.5)
net.w1<-cbind(w,b)
W<-2*matrix(runif(net.nhidden*net.nout)-0.5,net.nout,net.nhidden)
B<-2*(runif(net.nout)-0.5)
net.w2<-cbind(W,B)
traind_s<-traind
traind_s[,1]<-traind[,1]-mean(traind[,1])
traind_s[,2]<-traind[,2]-mean(traind[,2])
traind_s[,1]<-traind_s[,1]/sd(traind_s[,1])
traind_s[,2]<-traind_s[,2]/sd(traind_s[,2])
sampinex<-rbind(t(traind_s),rep(1,ntrainnum))
expectedout<-train1
eps<-0.01
a<-0.3
mc<-0.8
maxiter<-2000
iter<-0
errrec<-rep(0,maxiter)
outrec<-matrix(rep(0,ntrainnum*maxiter),ntrainnum,maxiter)
sigmoid<-function(x){
y<-1/(1+exp(-x))
return(y)
}
for(i in 1:maxiter){
hid_input<-net.w1%*%sampinex;
hid_out<-sigmoid(hid_input);
out_input1<-rbind(hid_out,rep(1,ntrainnum));
out_input2<-net.w2%*%out_input1;
out_out<-sigmoid(out_input2);
outrec[,i]<-t(out_out);
err<-expectedout-out_out;
sse<-sum(err^2);
errrec[i]<-sse;
iter<-iter+1;
if(sse<=eps)
break
Delta<-err*sigmoid(out_out)*(1-sigmoid(out_out))
delta<-(matrix(net.w2[,1:(length(net.w2[1,])-1)]))%*%Delta*sigmoid(hid_out)*(1-sigmoid(hid_out));
dWex<-Delta%*%t(out_input1)
dwex<-delta%*%t(sampinex)
if(i==1){
net.w2<-net.w2+a*dWex;
net.w1<-net.w1+a*dwex;
}
else{
net.w2<-net.w2+(1-mc)*a*dWex+mc*dWexold;
net.w1<-net.w1+(1-mc)*a*dwex+mc*dwexold;
}
dWexold<-dWex;
dwexold<-dwex;
}
testd_s<-testd
testd_s[,1]<-testd[,1]-mean(testd[,1])
testd_s[,2]<-testd[,2]-mean(testd[,2])
testd_s[,1]<-testd_s[,1]/sd(testd_s[,1])
testd_s[,2]<-testd_s[,2]/sd(testd_s[,2])
inex<-rbind(t(testd_s),rep(1,150-ntrainnum))
hid_input<-net.w1%*%inex
hid_out<-sigmoid(hid_input)
out_input1<-rbind(hid_out,rep(1,150-ntrainnum))
out_input2<-net.w2%*%out_input1
out_out<-sigmoid(out_input2)
out_out1<-out_out
out_out1[out_out<0.5]<-0
out_out1[out_out>=0.5]<-1
rate<-sum(out_out1==test1)/length(test1)
分类正确率为：0.9333333，是一个不错的学习器。这里须要注意的是动量因子mc的选定，mc不能过小，否则容易深陷局部最小而出不去人工神经网络和支持向量机优点，在本例中，如果mc=0.5，分类正确率仅为：0.5333333，学习疗效很不理想。
四、R中的神经网路函数
单层的前向神经网路模型在包nnet中的nnet函数，其调用格式为：
nnet(formula,data, weights, size, Wts, linout = F, entropy = F,
softmax = F, skip = F, rang = 0.7,decay = 0, maxit = 100,
trace = T)
参数说明:
size, 隐层结点数；
decay, 表明残差是递减的（可以避免过拟合）；
linout, 线性输出单元开关；
skip，是否准许跳过隐层；
maxit, 最大迭代次数；
Hess, 是否输出Hessian值
适用于神经网路的方式有predict,print和summary等，nnetHess函数拿来估算在考虑了权重参数下的Hessian矩阵，并且检验是否是局部最小。
我们使用nnet函数剖析Vehicle数据。随机选择半数观测作为训练集，剩下的作为测试集，构建只有包含3个节点的一个隐藏层的神经网路。输入如下程序：
library(nnet); #安装nnet软件包
library(mlbench); #安装mlbench软件包
data(Vehicle); #调入数据
n=length(Vehicle[,1]); #样本量
set.seed(1); #设随机数种子
samp=sample(1:n,n/2); #随机选择半数观测作为训练集
b=class.ind(Vehicle$Class); #生成类别的示性函数
test.cl=function(true,pred){true<-max.col(true);cres=max.col(pred);table(true,cres)};
a=nnet(Vehicle[samp,-19],b[samp,],size=3,rang=0.1,decay=5e-4,maxit=200); #利用训练集中前18个变量作为输入变量，隐藏层有3个节点，初始随机权值在[-0.1,0.1]，权值是逐渐衰减的。
test.cl(b[samp,],predict(a,Vehicle[samp,-19]))#给出训练集分类结果
test.cl(b[-samp,],predict(a,Vehicle[-samp,-19]));#给出测试集分类结果
#构建隐藏层包含15个节点的网络。接着上面的语句输入如下程序：
a=nnet(Vehicle[samp,-19],b[samp,],size=15,rang=0.1,decay=5e-4,maxit=10000);
test.cl(b[samp,],predict(a,Vehicle[samp,-19]));
test.cl(b[-samp,],predict(a,Vehicle[-samp,-19]));再看手写数字案例
最后，我们回到最开始的那种手写数字的案例，我们试着借助支持向量机重做这个案例。（这个案例的描述与数据参见《R语言与机器学习学习笔记（分类算法）（1）》）
由于nnet包对输入的维数有一定限制（我也不知道为何，可能在残差估算的时侯出现了一些bug，反正将支持向量机那一节的代码平行的移过来是会报错的）。我们这儿采用手写数字辨识技术中常用的办法处理这个案例：计算数字的特点。选择数字特点的办法有许多种，你随意百度一篇论文都有表述。我们这儿采用结构特点与统计特点结合的办法估算图象的特点。
我们这儿采用的统计特点与上图有一点的不同（结构特点一致），我们是将图片分为16块（4*4），统计每位方块中点的个数，这样我们就有25维的特点向量了。为了保证结果的可比性，我们也报告支持向量机的分类结果。
运行下述代码：
setwd("D:/R/data/digits/trainingDigits")
names<-list.files("D:/R/data/digits/trainingDigits")
data<-paste("train",1:1934,sep="")
for(i in 1:length(names))
assign(data[i],as.matrix(read.fwf(names[i],widths=rep(1,32))))
library(nnet)
label<-factor(rep(0:9,c(189,198,195,199,186,187,195,201,180,204)))
feature<-matrix(rep(0,length(names)*25),length(names),25)
for(i in 1:length(names)){
feature[i,1]<-sum(get(data[i])[,16])
feature[i,2]<-sum(get(data[i])[,8])
feature[i,3]<-sum(get(data[i])[,24])
feature[i,4]<-sum(get(data[i])[16,])
feature[i,5]<-sum(get(data[i])[11,])
feature[i,6]<-sum(get(data[i])[21,])
feature[i,7]<-sum(diag(get(data[i])))
feature[i,8]<-sum(diag(get(data[i])[,32:1]))
feature[i,9]<-sum((get(data[i])[17:32,17:32]))
feature[i,10]<-sum((get(data[i])[1:8,1:8]))
feature[i,11]<-sum((get(data[i])[9:16,1:8]))
feature[i,12]<-sum((get(data[i])[17:24,1:8]))
feature[i,13]<-sum((get(data[i])[25:32,1:8]))
feature[i,14]<-sum((get(data[i])[1:8,9:16]))
feature[i,15]<-sum((get(data[i])[9:16,9:16]))
feature[i,16]<-sum((get(data[i])[17:24,9:16]))
feature[i,17]<-sum((get(data[i])[25:32,9:16]))
feature[i,18]<-sum((get(data[i])[1:8,17:24]))
feature[i,19]<-sum((get(data[i])[9:16,17:24]))
feature[i,20]<-sum((get(data[i])[17:24,17:24]))
feature[i,21]<-sum((get(data[i])[25:32,17:24]))
feature[i,22]<-sum((get(data[i])[1:8,25:32]))
feature[i,23]<-sum((get(data[i])[9:16,25:32]))
feature[i,24]<-sum((get(data[i])[17:24,25:32]))
feature[i,25]<-sum((get(data[i])[25:32,25:32]))
}
data1 <- data.frame(feature,label)
m1<-nnet(label~.,data=data1,size=25,maxit = 2000,decay = 5e-6, rang = 0.1)
pred<-predict(m1,data1,type="class")
table(pred,label)
sum(diag(table(pred,label)))/length(names)
library("e1071")
m <- svm(feature,label,cross=10,type="C-classification")
m
summary(m)
pred<-fitted(m)
table(pred,label)
setwd("D:/R/data/digits/testDigits")
name<-list.files("D:/R/data/digits/testDigits")
data1<-paste("train",1:1934,sep="")
for(i in 1:length(name))
assign(data1[i],as.matrix(read.fwf(name[i],widths=rep(1,32))))
feature<-matrix(rep(0,length(name)*25),length(name),25)
for(i in 1:length(name)){
feature[i,1]<-sum(get(data1[i])[,16])
feature[i,2]<-sum(get(data1[i])[,8])
feature[i,3]<-sum(get(data1[i])[,24])
feature[i,4]<-sum(get(data1[i])[16,])
feature[i,5]<-sum(get(data1[i])[11,])
feature[i,6]<-sum(get(data1[i])[21,])
feature[i,7]<-sum(diag(get(data1[i])))
feature[i,8]<-sum(diag(get(data1[i])[,32:1]))
feature[i,9]<-sum((get(data1[i])[17:32,17:32]))
feature[i,10]<-sum((get(data1[i])[1:8,1:8]))
feature[i,11]<-sum((get(data1[i])[9:16,1:8]))
feature[i,12]<-sum((get(data1[i])[17:24,1:8]))
feature[i,13]<-sum((get(data1[i])[25:32,1:8]))
feature[i,14]<-sum((get(data1[i])[1:8,9:16]))
feature[i,15]<-sum((get(data1[i])[9:16,9:16]))
feature[i,16]<-sum((get(data1[i])[17:24,9:16]))
feature[i,17]<-sum((get(data1[i])[25:32,9:16]))
feature[i,18]<-sum((get(data1[i])[1:8,17:24]))
feature[i,19]<-sum((get(data1[i])[9:16,17:24]))
feature[i,20]<-sum((get(data1[i])[17:24,17:24]))
feature[i,21]<-sum((get(data1[i])[25:32,17:24]))
feature[i,22]<-sum((get(data1[i])[1:8,25:32]))
feature[i,23]<-sum((get(data1[i])[9:16,25:32]))
feature[i,24]<-sum((get(data1[i])[17:24,25:32]))
feature[i,25]<-sum((get(data1[i])[25:32,25:32]))
}
labeltest<-factor(rep(0:9,c(87,97,92,85,114,108,87,96,91,89)))
data2<-data.frame(feature,labeltest)
pred1<-predict(m1,data2,type="class")
table(pred1,labeltest)
sum(diag(table(pred1,labeltest)))/length(name)
pred<-predict(m,feature)
table(pred,labeltest)
sum(diag(table(pred,labeltest)))/length(name)
经整理，我们有如下输出结果：
可以看见，神经网路与支持向量机还是有一定的可比性，但支持向量机的结果还是要优于神经网路的。
这里我们神经网路取25个节点（隐藏层）似乎出现了过拟合的现象（虽然还不算过分严重）我们应当减轻节点个数得到更佳的预测结果。
关于节点的选择是个经验活，我们没有一定的规则。可以多试几次，结合训练集正确率与测试集正确率综合研判，但是构造神经网路的代价是昂贵的，所以有一个不太坏的结果也就可以停止了。（其他参数的选择同样这么，但是不如size那么重要）
特征的选定对于辨识问题来说相当的重要，也许主成份在选择特点时作用会比我们这样的选择更好，但是代价也更高，还有我们应当怎样选择主成份，怎么选择（选择哪张图的主成份）都是须要考虑的。
五、神经网路还是支持向量机
从里面的表述可以看出，神经网路与我们上面说的支持向量机有不少相像的地方，那么我们应当选择谁呢？下面是两种方式的一个简明对比：
正如费曼强调的那样“我们必须从一开始就澄清一个观点，就是假如某事不是科学，它并不一定不好。比如说，爱情就不是科学。因此，如果我们说某事不是科学，并不是说它有哪些不对，而只是说它不是科学。”与SVM相比，ANN不像一门科学，更像一门工程方法，但并不意味着它就一定就不好。
Furtherreading：
