R语言实现神经网络

人工神经网络（ANN），简称神经网络，是一种模仿生物神经网络的结构和功能的数学模型或计算模型。神经网络由大量的人工神经元联结进行计算。大多数情况下人工神经网络能在外界信息的基础上改变内部结构，是一种自适应系统。现代神经网络是一种非线性统计性数据建模工具，常用来对输入和输出间复杂的关系进行建模，或用来探索数据的模式。

     人工神经网络从以下四个方面去模拟人的智能行为:

• 物理结构：人工神经元将模拟生物神经元的功能

• 计算模拟：人脑的神经元有局部计算和存储的功能，通过连接构成一个系统。人工神经网络中也有大量有局部处理能力的神经元，也能够将信息进行大规模并行处理

• 存储与操作：人脑和人工神经网络都是通过神经元的连接强度来实现记忆存储功能，同时为概括、类比、推广提供有力的支持

• 训练：同人脑一样，人工神经网络将根据自己的结构特性，使用不同的训练、学习过程，自动从实践中获得相关知识

       神经网络是一种运算模型，由大量的节点（或称“神经元”，或“单元”）和之间相互联接构成。每个节点代表一种特定的输出函数，称为激励函数。每两个节点间的连接都代表一个对于通过该连接信号的加权值，称之为权重，这相当于人工神经网络的记忆。网络的输出则依网络的连接方式，权重值和激励函数的不同而不同。而网络自身通常都是对自然界某种算法或者函数的逼近，也可能是对一种逻辑策略的表达。

一、感知器

        感知器相当于神经网络的一个单层，由一个线性组合器和一个二值阈值原件构成：

     构成ANN系统的单层感知器：

• 感知器以一个实数值向量作为输入，计算这些输入的线性组合，如果结果大于某个阈值，就输出1，否则输出‐1。

• 感知器函数可写为：sign（w*x）有时可加入偏置b，写为sign（w*x b）

• 学习一个感知器意味着选择权w0,…,wn的值。所以感知器学习要考虑的候选假设空间H就是所有可能的实数值权向量的集合

     算法训练步骤：

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<a<1

6、判断，若满足收敛条件，算法结束，否则返回3

      注意，其中学习率a为了权值的稳定性不应过大，为了体现误差对权值的修正不应过小，说到底，这是个经验问题。

      从前面的叙述来看，感知器对于线性可分的例子是一定收敛的，对于不可分问题，它没法实现正确分类。这里与我们前面讲到的支持向量机的想法十分的相近，只是确定分类直线的办法有所不同。可以这么说，对于线性可分的例子，支持向量机找到了“最优的”那条分类直线，而单层感知器找到了一条可行的直线。

      我们以鸢尾花数据集为例，由于单层感知器是一个二分类器，所以我们将鸢尾花数据也分为两类，“setosa”与“versicolor”（将后两类均看做第2类），那么数据按照特征：花瓣长度与宽度做分类。

运行下面的代码：

[plain] view plain copy print ?

1. #感知器训练结果：  

2. a<-0.2  

3. w<-rep(0,3)  

4. iris1<-t(as.matrix(iris[,3:4]))  

5. d<-c(rep(0,50),rep(1,100))  

6. e<-rep(0,150)  

7. p<-rbind(rep(1,150),iris1)  

8. max<-100000  

9. eps<-rep(0,100000)  

10. i<-0  

11. repeat{  

12.    v<-w%*%p;  

13.    y<-ifelse(sign(v)>=0,1,0);  

14.    e<-d-y;  

15.    eps[i 1]<-sum(abs(e))/length(e)  

16.    if(eps[i 1]<0.01){  

17.        print("finish:");  

18.        print(w);  

19.        break;  

20.    }  

21.    w<-w a*(d-y)%*%t(p);  

22.    i<-i 1;  

23.    if(i>max){  

24.        print("max time loop");  

25.        print(eps[i])  

26.        print(y);  

27.        break;  

28.    }  

29. }  

30. #绘图程序  

31. plot(Petal.Length~Petal.Width,xlim=c(0,3),ylim=c(0,8),  

32. data=iris[iris$Species=="virginica",])  

33. data1<-iris[iris$Species=="versicolor",]  

34. points(data1$Petal.Width,data1$Petal.Length,col=2)  

35. data2<-iris[iris$Species=="setosa",]  

36. points(data2$Petal.Width,data2$Petal.Length,col=3)  

37. x<-seq(0,3,0.01)  

38. y<-x*(-w[2]/w[3])-w[1]/w[3]  

39. lines(x,y,col=4)  

40. #绘制每次迭代的平均绝对误差  

41. plot(1:i,eps[1:i],type="o")  

      分类结果如图：

       这是运行了7次得到的结果。与我们前面的支持向量机相比，显然神经网络的单层感知器分类不是那么的可信，有些弱。

      我们可以尝试来做交叉验证，可以发现交叉验证结果并不理想。

二、线性神经网络

      尽管当训练样例线性可分时，感知器法则可以成功地找到一个权向量，但如果样例不是线性可分时它将不能收敛。因此，人们设计了另一个训练法则来克服这个不足，称为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<2/maxeign, delta法则按方差收敛。其中maxeigen为x’x的最大特征值。故我们这里使用1/maxeign作为a的值。

       我们还是以上面的鸢尾花数据为例来说这个问题。运行代码：

[plain] view plain copy print ?

1. p<-rbind(rep(1,150),iris1)  

2. d<-c(rep(0,50),rep(1,100))  

3. w<-rep(0,3)  

4. a<-1/max(eigen(t(p)%*%p)$values)  

5. max<-1000  

6. e<-rep(0,150)  

7. eps<-rep(0,1000)  

8. i<-0  

9. for(i in 1:max){  

10.    v<-w%*%p;  

11.    y<-v;  

12.    e<-d-y;  

13.    eps[i 1]<-sum(e^2)/length(e)  

14.    w<-w a*(d-y)%*%t(p);  

15.    if(i==max)  

16.        print(w)  

17. }  

       得到分类直线：

      相比感知器分类而言已经好了太多了，究其原因不外乎传递函数由二值阈值函数变为了线性函数，这也就是我们前面提到的delta法则会收敛到目标概念的最佳近似。增量法则渐近收敛到最小误差假设，可能需要无限的时间，但无论训练样例是否线性可分都会收敛。

       为了明了这一点我们考虑鸢尾花数据后两类花的分类（这里我们将前两类看做一类），使用感知器：

使用线性分类器：

       但是要解释的一点是，收敛并不意味着分类效果更好，要解决线性不可分问题需要的是添加非线性输入或者增加神经元。我们以Minsky & Papert (1969)提出的异或例子为例说明这一点。

使用线性神经网络，代码与上面完全相同，略。

第一个神经元输出：

权值： [,1]     [,2]    [,3]

[1,]    0.75   0.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  TRUE TRUE   FALSE

即0，1，1，0，分类正确。

       最后再说一点，Delta规则只能训练单层网络，但这不会对其功能造成很大的影响。从理论上说，多层神经网络并不比单层神经网络更强大，他们具有同样的能力。

三、BP神经网络

1、sigmoid函数分类

        回顾我们前面提到的感知器，它使用示性函数作为分类的办法。然而示性函数作为分类器它的跳点让人觉得很难处理，幸好sigmoid函数y=1/(1 e^-x)有类似的性质，且有着光滑性这一优良性质。我们通过下图可以看见sigmoid函数的图像：

Sigmoid函数有着计算代价不高，易于理解与实现的优点但也有着欠拟合，分类精度不高的特性，我们在支持向量机一章中就可以看到sigmoid函数差劲的分类结果。

2、BP神经网络结构

       BP (Back Propagation)神经网络，即误差反传误差反向传播算法的学习过程，由信息的正向传播和误差的反向传播两个过程组成。由下图可知，BP神经网络是一个三层的网络:

• 输入层(input layer)：输入层各神经元负责接收来自外界的输入信息，并传递给中间层各神经元；

• 隐藏层(Hidden Layer):中间层是内部信息处理层，负责信息变换，根据信息变化能力的需求，中间层可以设计为单隐层或者多隐层结构；最后一个隐层传递到输出层各神经元的信息，经进一步处理后，完成一次学习的正向传播处理过程；

• 输出层(Output Layer):顾名思义，输出层向外界输出信息处理结果；

       当实际输出与期望输出不符时，进入误差的反向传播阶段。误差通过输出层，按误差梯度下降的方式修正各层权值，向隐藏层、输入层逐层反传。周而复始的信息正向传播和误差反向传播过程，是各层权值不断调整的过程，也是神经网络学习训练的过程，此过程一直进行到网络输出的误差减少到可以接受的程度，或者预先设定的学习次数为止。

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 造成网络整体的误差也越多。另一方面，如果输出节点影响网络整体的误差越少，n 对输出节点的影响也相应减少。这里d(j)是对网络的整体误差的基值，W(n,j) 是 n 对 j 造成的影响，d(j) * W(n,j) 是这两种影响的总和。但是 n 几乎总是影响多个输出节点，也许会影响每一个输出结点，这样，d(n) 可以表示为：SUM(d(j)*W(n,j))

       这里j是一个从n获得输入的输出节点，联系起来，我们就得到了一个培训规则。

• 第1部分：在隐藏节点n和输出节点j之间权系数改变，如下所示：
  alpha *s'(a(p,n))*(t(p,n) - y(p,n)) * X(p,n,j)

• 第 2 部分：在输入节点i和输出节点n之间权系数改变，如下所示：
  alpha *s'(a(p,n)) * sum(d(j) * W(n,j)) * X(p,i,n)

      这里每个从n接收输入的输出节点j都不同。关于反向传播算法的基本情况大致如此。

      通常把第 1部分称为正向传播，把第2部分称为反向传播。反向传播的名字由此而来。

4、最速下降法与其改进

      最速下降法的基本思想是：要找到某函数的最小值，最好的办法是沿函数的梯度方向探寻，如果梯度记为d,那么迭代公式可写为w=w-alpha*d，其中alpha可理解为我们前面提到的学习速率。

      最速下降法有着收敛速度慢（因为每次搜索与前一次均正交，收敛是锯齿形的），容易陷入局部最小值等缺点，所以他的改进办法也有不少，最常见的是增加动量项与学习率可变。

     增加冲量项（Momentum）

• 修改权值更新法则，使第n次迭代时的权值的更新部分地依赖于发生在第n‐1次迭代时的更新
  Delta（w）(n)=-alpha*(1-mc)*Delta(w)(n) mc*Delta(w)(n-1)

• 右侧第一项就是权值更新法则，第二项被称为冲量项

• 梯度下降的搜索轨迹就像一个球沿误差曲面滚下，冲量使球从一次迭代到下一次迭代时以同样的方向滚动

• 冲量有时会使这个球滚过误差曲面的局部极小值或平坦区域

• 冲量也具有在梯度不变的区域逐渐增大搜索步长的效果，从而加快收敛。

     改变学习率

• 当误差减小趋近目标时，说明修正方向是正确的，可以增加学习率；当误差增加超过一个范围时，说明修改不正确，需要降低学习率。

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代码：

[plain] view plain copy print ?

1. iris1<-as.matrix(iris[,3:4])  

2. iris1<-cbind(iris1,c(rep(1,100),rep(0,50)))  

3. set.seed(5)  

4. n<-length(iris1[,1])  

5. samp<-sample(1:n,n/5)  

6. traind<-iris1[-samp,c(1,2)]  

7. train1<-iris1[-samp,3]  

8. testd<-iris1[samp,c(1,2)]  

9. test1<-iris1[samp,3]  

10.  

11. set.seed(1)  

12. ntrainnum<-120  

13. nsampdim<-2  

14.  

15. net.nin<-2  

16. net.nhidden<-3  

17. net.nout<-1  

18. w<-2*matrix(runif(net.nhidden*net.nin)-0.5,net.nhidden,net.nin)  

19. b<-2*(runif(net.nhidden)-0.5)  

20. net.w1<-cbind(w,b)  

21. W<-2*matrix(runif(net.nhidden*net.nout)-0.5,net.nout,net.nhidden)  

22. B<-2*(runif(net.nout)-0.5)  

23. net.w2<-cbind(W,B)  

24.  

25. traind_s<-traind  

26. traind_s[,1]<-traind[,1]-mean(traind[,1])  

27. traind_s[,2]<-traind[,2]-mean(traind[,2])  

28. traind_s[,1]<-traind_s[,1]/sd(traind_s[,1])  

29. traind_s[,2]<-traind_s[,2]/sd(traind_s[,2])  

30.  

31. sampinex<-rbind(t(traind_s),rep(1,ntrainnum))  

32. expectedout<-train1  

33.  

34. eps<-0.01  

35. a<-0.3  

36. mc<-0.8  

37. maxiter<-2000  

38. iter<-0  

39.  

40. errrec<-rep(0,maxiter)  

41. outrec<-matrix(rep(0,ntrainnum*maxiter),ntrainnum,maxiter)  

42.  

43. sigmoid<-function(x){  

44.    y<-1/(1 exp(-x))  

45.    return(y)  

46. }  

47.  

48. for(i in 1:maxiter){  

49.    hid_input<-net.w1%*%sampinex;  

50.    hid_out<-sigmoid(hid_input);  

51.    out_input1<-rbind(hid_out,rep(1,ntrainnum));  

52.    out_input2<-net.w2%*%out_input1;  

53.    out_out<-sigmoid(out_input2);  

54.    outrec[,i]<-t(out_out);  

55.    err<-expectedout-out_out;  

56.    sse<-sum(err^2);  

57.    errrec[i]<-sse;  

58.    iter<-iter 1;  

59.    if(sse<=eps)  

60.        break  

61.      

62.    Delta<-err*sigmoid(out_out)*(1-sigmoid(out_out))  

63.    delta<-(matrix(net.w2[,1:(length(net.w2[1,])-1)]))%*%Delta*sigmoid(hid_out)*(1-sigmoid(hid_out));  

64.  

65.    dWex<-Delta%*%t(out_input1)  

66.    dwex<-delta%*%t(sampinex)  

67.  

68.    if(i==1){  

69.            net.w2<-net.w2 a*dWex;  

70.            net.w1<-net.w1 a*dwex;  

71.        }  

72.    else{  

73.            net.w2<-net.w2 (1-mc)*a*dWex mc*dWexold;  

74.            net.w1<-net.w1 (1-mc)*a*dwex mc*dwexold;  

75.    }  

76.  

77.    dWexold<-dWex;  

78.    dwexold<-dwex;  

79. }  

80.  

81.  

82. testd_s<-testd  

83. testd_s[,1]<-testd[,1]-mean(testd[,1])  

84. testd_s[,2]<-testd[,2]-mean(testd[,2])  

85. testd_s[,1]<-testd_s[,1]/sd(testd_s[,1])  

86. testd_s[,2]<-testd_s[,2]/sd(testd_s[,2])  

87.  

88. inex<-rbind(t(testd_s),rep(1,150-ntrainnum))  

89. hid_input<-net.w1%*%inex  

90. hid_out<-sigmoid(hid_input)  

91. out_input1<-rbind(hid_out,rep(1,150-ntrainnum))  

92. out_input2<-net.w2%*%out_input1  

93. out_out<-sigmoid(out_input2)  

94. out_out1<-out_out  

95.  

96. out_out1[out_out<0.5]<-0  

97. out_out1[out_out>=0.5]<-1  

98.  

99. 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个节点的一个隐藏层的神经网络。输入如下程序：

[plain] view plain copy print ?

1. library(nnet);  #安装nnet软件包  

2. library(mlbench);  #安装mlbench软件包  

3. data(Vehicle);  #调入数据  

4. n=length(Vehicle[,1]); #样本量  

5. set.seed(1);  #设随机数种子  

6. samp=sample(1:n,n/2);  #随机选择半数观测作为训练集  

7. b=class.ind(Vehicle$Class);  #生成类别的示性函数  

8. test.cl=function(true,pred){true<-max.col(true);cres=max.col(pred);table(true,cres)};  

9. a=nnet(Vehicle[samp,-19],b[samp,],size=3,rang=0.1,decay=5e-4,maxit=200);  #利用训练集中前18个变量作为输入变量，隐藏层有3个节点，初始随机权值在[-0.1,0.1]，权值是逐渐衰减的。  

10. test.cl(b[samp,],predict(a,Vehicle[samp,-19]))#给出训练集分类结果  

11. test.cl(b[-samp,],predict(a,Vehicle[-samp,-19]));#给出测试集分类结果  

12. #构建隐藏层包含15个节点的网络。接着上面的语句输入如下程序：  

13. a=nnet(Vehicle[samp,-19],b[samp,],size=15,rang=0.1,decay=5e-4,maxit=10000);  

14. test.cl(b[samp,],predict(a,Vehicle[samp,-19]));  

15. test.cl(b[-samp,],predict(a,Vehicle[-samp,-19]));  

再看手写数字案例

       最后，我们回到最开始的那个手写数字的案例，我们试着利用支持向量机重做这个案例。（这个案例的描述与数据参见《R语言与机器学习学习笔记（分类算法）（1）》）

      由于nnet包对输入的维数有一定限制（我也不知道为什么，可能在权值计算的时候出现了一些bug，反正将支持向量机那一节的代码平行的移过来是会报错的）。我们这里采用手写数字识别技术中常用的办法处理这个案例：计算数字的特征。选择数字特征的办法有许多种，你随便百度一篇论文都有叙述。我们这里采用结构特征与统计特征结合的办法计算图像的特征。

     我们这里采用的统计特征与上图有一点的不同（结构特征一致），我们是将图片分为16块（4*4），统计每个小方块中点的个数，这样我们就有25维的特征向量了。为了保证结果的可比性，我们也报告支持向量机的分类结果。

      运行下列代码：

[plain] view plain copy print ?

1. setwd("D:/R/data/digits/trainingDigits")  

2. names<-list.files("D:/R/data/digits/trainingDigits")  

3. data<-paste("train",1:1934,sep="")  

4. for(i in 1:length(names))  

5.         assign(data[i],as.matrix(read.fwf(names[i],widths=rep(1,32))))  

6. library(nnet)  

7. label<-factor(rep(0:9,c(189,198,195,199,186,187,195,201,180,204)))  

8.  

9. feature<-matrix(rep(0,length(names)*25),length(names),25)  

10. for(i in 1:length(names)){  

11.        feature[i,1]<-sum(get(data[i])[,16])  

12.        feature[i,2]<-sum(get(data[i])[,8])  

13.        feature[i,3]<-sum(get(data[i])[,24])  

14.        feature[i,4]<-sum(get(data[i])[16,])  

15.        feature[i,5]<-sum(get(data[i])[11,])  

16.        feature[i,6]<-sum(get(data[i])[21,])  

17.        feature[i,7]<-sum(diag(get(data[i])))  

18.        feature[i,8]<-sum(diag(get(data[i])[,32:1]))  

19.        feature[i,9]<-sum((get(data[i])[17:32,17:32]))  

20.        feature[i,10]<-sum((get(data[i])[1:8,1:8]))  

21.        feature[i,11]<-sum((get(data[i])[9:16,1:8]))  

22.        feature[i,12]<-sum((get(data[i])[17:24,1:8]))  

23.        feature[i,13]<-sum((get(data[i])[25:32,1:8]))  

24.        feature[i,14]<-sum((get(data[i])[1:8,9:16]))  

25.        feature[i,15]<-sum((get(data[i])[9:16,9:16]))  

26.        feature[i,16]<-sum((get(data[i])[17:24,9:16]))  

27.        feature[i,17]<-sum((get(data[i])[25:32,9:16]))  

28.        feature[i,18]<-sum((get(data[i])[1:8,17:24]))  

29.        feature[i,19]<-sum((get(data[i])[9:16,17:24]))  

30.        feature[i,20]<-sum((get(data[i])[17:24,17:24]))  

31.        feature[i,21]<-sum((get(data[i])[25:32,17:24]))  

32.        feature[i,22]<-sum((get(data[i])[1:8,25:32]))  

33.        feature[i,23]<-sum((get(data[i])[9:16,25:32]))  

34.        feature[i,24]<-sum((get(data[i])[17:24,25:32]))  

35.        feature[i,25]<-sum((get(data[i])[25:32,25:32]))  

36. }  

37. data1 <- data.frame(feature,label)  

38. m1<-nnet(label~.,data=data1,size=25,maxit = 2000,decay = 5e-6, rang = 0.1)  

39. pred<-predict(m1,data1,type="class")  

40. table(pred,label)  

41. sum(diag(table(pred,label)))/length(names)  

42.  

43. library("e1071")  

44. m <- svm(feature,label,cross=10,type="C-classification")  

45. m  

46. summary(m)  

47. pred<-fitted(m)  

48. table(pred,label)  

49.  

50. setwd("D:/R/data/digits/testDigits")  

51. name<-list.files("D:/R/data/digits/testDigits")  

52. data1<-paste("train",1:1934,sep="")  

53. for(i in 1:length(name))  

54.    assign(data1[i],as.matrix(read.fwf(name[i],widths=rep(1,32))))  

55.  

56. feature<-matrix(rep(0,length(name)*25),length(name),25)  

57. for(i in 1:length(name)){  

58.        feature[i,1]<-sum(get(data1[i])[,16])  

59.        feature[i,2]<-sum(get(data1[i])[,8])  

60.        feature[i,3]<-sum(get(data1[i])[,24])  

61.        feature[i,4]<-sum(get(data1[i])[16,])  

62.        feature[i,5]<-sum(get(data1[i])[11,])  

63.        feature[i,6]<-sum(get(data1[i])[21,])  

64.        feature[i,7]<-sum(diag(get(data1[i])))  

65.        feature[i,8]<-sum(diag(get(data1[i])[,32:1]))  

66.        feature[i,9]<-sum((get(data1[i])[17:32,17:32]))  

67.        feature[i,10]<-sum((get(data1[i])[1:8,1:8]))  

68.        feature[i,11]<-sum((get(data1[i])[9:16,1:8]))  

69.        feature[i,12]<-sum((get(data1[i])[17:24,1:8]))  

70.        feature[i,13]<-sum((get(data1[i])[25:32,1:8]))  

71.        feature[i,14]<-sum((get(data1[i])[1:8,9:16]))  

72.        feature[i,15]<-sum((get(data1[i])[9:16,9:16]))  

73.        feature[i,16]<-sum((get(data1[i])[17:24,9:16]))  

74.        feature[i,17]<-sum((get(data1[i])[25:32,9:16]))  

75.        feature[i,18]<-sum((get(data1[i])[1:8,17:24]))  

76.        feature[i,19]<-sum((get(data1[i])[9:16,17:24]))  

77.        feature[i,20]<-sum((get(data1[i])[17:24,17:24]))  

78.        feature[i,21]<-sum((get(data1[i])[25:32,17:24]))  

79.        feature[i,22]<-sum((get(data1[i])[1:8,25:32]))  

80.        feature[i,23]<-sum((get(data1[i])[9:16,25:32]))  

81.        feature[i,24]<-sum((get(data1[i])[17:24,25:32]))  

82.        feature[i,25]<-sum((get(data1[i])[25:32,25:32]))  

83. }  

84. labeltest<-factor(rep(0:9,c(87,97,92,85,114,108,87,96,91,89)))  

85. data2<-data.frame(feature,labeltest)  

86. pred1<-predict(m1,data2,type="class")  

87. table(pred1,labeltest)  

88. sum(diag(table(pred1,labeltest)))/length(name)  

89.  

90. pred<-predict(m,feature)  

91. table(pred,labeltest)  

92. sum(diag(table(pred,labeltest)))/length(name)  

      经整理，我们有如下输出结果：

       可以看到，神经网络与支持向量机还是有一定的可比性，但支持向量机的结果还是要优于神经网络的。

       这里我们神经网络取25个节点（隐藏层）似乎出现了过拟合的现象（虽然还不算过于严重）我们应该减少节点个数得到更佳的预测结果。

       关于节点的选择是个经验活，我们没有一定的规则。可以多试几次，结合训练集正确率与测试集正确率综合研判，但是构造神经网络的代价是高昂的，所以有一个不太坏的结果也就可以停止了。（其他参数的选择同样如此，但是不如size那么重要）    

      特征的选取对于识别问题来说相当的重要，也许主成分在选择特征时作用会比我们这样的选择更好，但是代价也更高，还有我们应该如何选择主成分，怎么选择（选择哪张图的主成分）都是需要考虑的。

五、神经网络还是支持向量机

      从上面的叙述可以看出，神经网络与我们前面说的支持向量机有不少相似的地方，那么我们应该选择谁呢？下面是两种方法的一个简明对比：

• – SVM的理论基础比NN更坚实，更像一门严谨的“科学”（三要素：问题的表示、问题的解决、证明）

• – SVM ——严格的数学推理
  –ANN ——强烈依赖于工程技巧

• –推广能力取决于“经验风险值”和“置信范围值”，ANN不能控制两者中的任何一个。

• –ANN设计者用高超的工程技巧弥补了数学上的缺陷——设计特殊的结构，利用启发式算法，有时能得到出人意料的好结果。

       正如费曼指出的那样“我们必须从一开始就澄清一个观点，就是如果某事不是科学，它并不一定不好。比如说，爱情就不是科学。因此，如果我们说某事不是科学，并不是说它有什么不对，而只是说它不是科学。”与SVM相比，ANN不像一门科学，更像一门工程技巧，但并不意味着它就一定就不