这个也是本学期模式识别的一个小作业题。
实验目的：采用神经网络的方法对二维双螺旋样本点进行分类。
实验数据：由于本次的实验数据是双螺旋，需要用数学公式产生，其产生方法见参考文献[1].
即由下面的公式产生：
实验数据分为2部分：训练数据和测试数据。由于双螺旋曲线有自己的方程表达式，我们产生的测试数据不能与训练数据一样，否则训练出来的网络没有说服力，因此我们在程序中2个样本集在双螺旋曲线上的采样间隔不同，且起始位置不同，这样保证了测试样本和训练样本的数据集没有重叠，训练出来的模型的分类能力更具有说服力。
下图是试验过程中产生的训练样本和测试样本分布图：
其中训练样本200个，每条螺旋线100个；测试样本数据也为200个，每天螺旋线100个。
试验方法：本次试验我采用了2种方法，第一种是利用matlab自带的ANN工具箱函数，来训练出一个BP神经网络；另一种方法是不采用matlab自带的工具箱，而自己用代码实现 一个BP神经网络，最后给出了2者的对比。
方法一：采用matlab自带神经网络工具箱
实验说明：采用的是matlab自带的BP神经网络，隐含层只有1层，该层有8个神经元。程序代码中的studytra_data, studytra_label为训练数据，分别是2*200的矩阵。 predict_test_data, predict_test_label 为测试数据，也是2*200的矩阵。网络形成用的是feedforwardnet函数，BP网络训练采用的是train函数。
实验结果：
实验预测的结果如下所示：
即第一条螺旋线的预测精度为92%,第二条螺旋线的预测精度为68%.
程序代码(matlab)：
1 %% 产生双螺旋数据 2 3 studytra_num=100; %训练样本数目 4 predict_test_number=100; %测试样本数目 5 6 i=(1:1:studytra_num)'; 7 %双螺旋数据点的产生方程 8 alpha1=pi*(i-1)/25; 9 beta=0.4*((105-i)/104); 10 x0=0.5+beta.*sin(alpha1); 11 y0=0.5+beta.*cos(alpha1); 12 z0=zeros(studytra_num,1); 13 x1=0.5-beta.*sin(alpha1); 14 y1=0.5-beta.*cos(alpha1); 15 z1=ones(studytra_num,1); 16 17 %% 事实证明BP神经网络在训练数据时与输入数据正负样本的顺序是有关系的 18 % 如果一开始的一半数据都是正(负)样本，后面的全是负(正)样本，则训练出来的 19 %效果不好，所以这里需要随机打乱 20 k=rand(1,2*studytra_num); 21 [m,n]=sort(k); 22 23 studytra=[x0 y0 z0;x1,y1,z1]; %1条螺旋线数据点,200*3的矩阵 24 trian_label1=studytra(n(1:2*studytra_num),end)'; %训练数据类别，1*200的行向量 25 studytra_data1=studytra(n(1:2*studytra_num),1:end-1)'; %训练数据属性，2*200的矩阵 26 27 %把1维的输出变成2维的输出,studytra_labe2为200*2的矩阵 28 for i=1:2*studytra_num 29 switch trian_label1(i) 30 case 0 31 studytra_label2(i,:)=[1 0]; 32 case 1 33 studytra_label2(i,:)=[0 1]; 34 end 35 end 36 37 [studytra_data,studytra_datas]=mapminmax(studytra_data1);%studytra_data为2*200的矩阵 38 studytra_label=studytra_label2'; %studytra_label为2*200的矩阵 39 40 plot(x0,y0,'r+'); 41 hold on; 42 plot(x1,y1,'go'); 43 44 45 %% 产生双螺旋测试数据 46 47 i=(1.5:1:predict_test_number+0.5)'; %每类51个样本 48 %双螺旋数据点的产生方程 49 alpha2=pi*(i-1)/25; 50 beta2=0.4*((105-i)/104); 51 m0=0.5+beta2.*sin(alpha2); 52 n0=0.5+beta2.*cos(alpha2); 53 s0=zeros(predict_test_number,1); 54 m1=0.5-beta2.*sin(alpha2); 55 n1=0.5-beta2.*cos(alpha2); 56 s1=ones(predict_test_number,1); 57 58 predict_test=[m0 n0 s0;m1,n1,s1]; %1条螺旋线数据点,3*102的矩阵 59 predict_test_label1=predict_test(:,end)'; %测试数据类别，1*102的行向量 60 predict_test_data1=predict_test(:,1:end-1)'; %测试数据属性，2*102的矩阵 61 62 %把1维的输出变成2维的输出,studytra_labe2为200*2的矩阵 63 for i=1:2*predict_test_number 64 switch predict_test_label1(i) 65 case 0 66 predict_test_label2(i,:)=[1 0]; 67 case 1 68 predict_test_label2(i,:)=[0 1]; 69 end 70 end 71 72 predict_test_label=predict_test_label2'; %predict_test_label为2*102的矩阵 73 74 %% 画出测试数据双螺旋曲线 75 plot(m0,n0,'c+'); 76 hold on; 77 plot(m1,n1,'yo'); 78 legend('训练数据螺旋线1','训练数据螺旋线2','测试数据螺旋线1','测试数据螺旋线2'); 79 80 predict_test_data=mapminmax('apply',predict_test_data1,studytra_datas); 81 82 %% 至此，训练和测试数据如下所示： 83 studytra_data ; %2*200的矩阵 84 studytra_label ; %2*200的矩阵 85 predict_test_data ; %2*200的矩阵 86 predict_test_label ; %2*200的矩阵 87 88 % %% 调用matlab的神经网络工具箱函数用来训练,方法1 89 % net=newff(minmax(studytra_data),[10,2],{'tansig','purelin','studytragdm'}) 90 % inputWeights=net.IW{1,1}; 91 % inputbias=net.b{1}; 92 % layerWeights=net.LW{2,1}; 93 % layerbias=net.b{2}; 94 % 95 % net.studytraParam.show=50; 96 % net.studytraParam.lr=0.05; 97 % net.studytraParam.mc=0.9; 98 % net.studytraParam.epochs=1000; 99 % net.studytraParam.goal=1e-2; 100 % 101 % [net,tr]=studytra(net,studytra_data,studytra_label); 102 103 104 %% 调用matlab的神经网络工具箱函数用来训练,方法2 105 net=feedforwardnet(8); 106 %注意此时进入train函数的样本每一列为一个样本的属性，列数为样本数，且命名时最好不要含有train,test字样，否则matlab会报错 107 net=train(net,studytra_data,studytra_label); 108 %view(net); 109 predict_label=sim(net,predict_test_data); 110 111 112 % %% 用训练到的模型预测数据 113 % for i=1:2*predict_test_number 114 % for j=1:midnum 115 % I(j)=predict_test_data(:,i)'*w1(j,:)'+b1(j); 116 % Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量 117 % end 118 % predict_test(:,i)=w2'*Iout'+b2;%predict_test为2*102的矩阵 119 % end 120 121 122 %% 预测结果分析 123 for i=1:2*predict_test_number 124 output_pred(i)=find(predict_label(:,i)==max(predict_label(:,i))); %out_pred为1*102的矩阵 125 end 126 127 error=output_pred-predict_test_label1-1; % 128 129 130 %% 计算出每一类预测错误的个数总和 131 k=zeros(1,2); %k=[0 0] 132 for i=1:2*predict_test_number 133 if error(i)~=0 %matlab中不能用if error(i)！=0 134 [b c]=max(predict_test_label(:,i)); 135 switch c 136 case 1 137 k(1)=k(1)+1; 138 case 2 139 k(2)=k(2)+1; 140 end 141 end 142 end 143 144 145 %% 求出每一类总体的个数和 146 kk=zeros(1,2); %k=[0 0] 147 for i=1:2*predict_test_number 148 [b c]=max(predict_test_label(:,i)); 149 switch c 150 case 1 151 kk(1)=kk(1)+1; 152 case 2 153 kk(2)=kk(2)+1; 154 end 155 end 156 157 158 %% 计算每一类的预测正确率 159 accuracy=(kk-k)./kk
方法二：自己实现BP神经网络的训练和测试函数
实验说明：实现BP神经网络的方法参考的是文献[2]，主要分为下面6个步骤：
- 网络的初始化
- 隐含层输出计算
- 输出层输出计算
- 误差计算
- 权值更新
- 阈值更新
同样，训练和测试样本都为2*200的矩阵，且训练过程中设置的循环次数为200次，只有1个隐含层，且神经元个数为5个。
实验结果：
实验预测的精度如下所示：
即第1条螺旋线的预测精度为94%,，第2条螺旋线的预测精度为37%。
实验代码(matlab)：
1 %% 产生双螺旋数据,每类100个样本点，共200个样本 2 train_num=100; 3 train_circle_number=5000; 4 test_number=100; 5 i=(1:1:train_num)'; 6 7 %双螺旋数据点的产生方程 8 alpha1=pi*(i-1)/25; 9 beta=0.4*((105-i)/104); 10 x0=0.5+beta.*sin(alpha1); 11 y0=0.5+beta.*cos(alpha1); 12 z0=zeros(train_num,1); 13 x1=0.5-beta.*sin(alpha1); 14 y1=0.5-beta.*cos(alpha1); 15 z1=ones(train_num,1); 16 17 %% 事实证明BP神经网络在训练数据时与输入数据正负样本的顺序是有关系的 18 % 如果一开始的一半数据都是正(负)样本，后面的全是负(正)样本，则训练出来的 19 %效果不好，所以这里需要随机打乱 20 k=rand(1,2*train_num); 21 [m,n]=sort(k); 22 23 train=[x0 y0 z0;x1,y1,z1]; %1条螺旋线数据点,200*3的矩阵 24 trian_label1=train(n(1:2*train_num),end)'; %训练数据类别，1*200的行向量 25 train_data1=train(n(1:2*train_num),1:end-1)'; %训练数据属性，2*200的矩阵 26 27 %把1维的输出变成2维的输出,train_labe2为200*2的矩阵 28 for i=1:2*train_num 29 switch trian_label1(i) 30 case 0 31 train_label2(i,:)=[1 0]; 32 case 1 33 train_label2(i,:)=[0 1]; 34 end 35 end 36 37 train_label=train_label2'; %train_label为2*200的矩阵 38 39 plot(x0,y0,'r+'); 40 hold on; 41 plot(x1,y1,'go'); 42 %legend(); 43 44 %% BP神经网络结构的初始化 45 %网络结构，2个输入，3个神经元，2个输出 46 innum=2; 47 midnum=5; 48 outnum=2; 49 50 [train_data,train_datas]=mapminmax(train_data1); 51 52 %输入输出取值阈值随机初始化 53 %w1矩阵表示每一行为一个隐含层神经元的输入权值 54 w1=rands(midnum,innum);%rands函数用来初始化神经元的权值和阈值是很合适的,w1为3*2的矩阵 55 b1=rands(midnum,1);%b1为3*1的矩阵 56 %w2矩阵表示每一列为一个输出层神经元的输入权值 57 w2=rands(midnum,outnum);%w2为3*2的矩阵 58 b2=rands(outnum,1);%b2为2*1的矩阵 59 60 %用来保存上一次的权值和阈值，因为后面的更新方差是递归的，要用到 61 w1_1=w1;w1_2=w1_1; 62 b1_1=b1;b1_2=b1_1; 63 w2_1=w2;w2_2=w2_1; 64 b2_1=b2;b2_2=b2_1; 65 66 %学习率的设定 67 alpha=0.05; 68 69 %训练10次就ok了，而不管训练后的结果如何 70 for train_circle=1:train_circle_number ; 71 for i=1:2*train_num; %200个训练样本 72 %% 输入层的输出 73 x=train_data(:,i);%取出第i个样本，x(i)为2*1的列向量 74 %% 隐含层的输出 75 for j=1:midnum; 76 I(j)=train_data(:,i)'*w1(j,:)'+b1(j); %I(j)为1*1的实数 77 Iout(j)=1/(1+exp(-I(j))); %Iout(j)也为1*1的实数 78 end %Iout为1*3的行向量 79 %% 输出层的输出 80 yn=(Iout*w2)'+b2; %yn为2*1的列向量，因此此时的传函为线性的，所以可以一步到位，不必上面 81 82 %% 计算误差 83 e=train_label(:,i)-yn; %e为2*1的列向量，保存的是误差值 84 85 %计算权值变换率 86 dw2=e*Iout; %dw2为2*3的矩阵，每一行表示输出接点的输入权值变化率 87 db2=e'; %e为1*2的行向量 88 89 for j=1:midnum 90 S=1/(1+exp(-I(j))); 91 FI(j)=S*(1-S); %FI(j)为一实数，FI为1*3的行向量 92 end 93 94 for k=1:1:innum 95 for j=1:midnum 96 dw1(k,j)=FI(j)*x(k)*(e(1)*w2(j,1)+e(2)*w2(j,2)); %dw1为2*3的矩阵 97 db1(j)=FI(j)*(e(1)*w2(j,1)+e(2)*w2(j,2)); %db1为1*3的矩阵 98 end 99 end 100 101 %% 权值更新方程 102 w1=w1_1+alpha*dw1'; %w1仍为3*2的矩阵 103 b1=b1_1+alpha*db1'; %b1仍为3*1的矩阵 104 w2=w2_1+alpha*dw2'; %w2仍为3*2的矩阵 105 b2=b2_1+alpha*db2'; %b2仍为2*1的矩阵 106 107 %% 保存上一次的权值和阈值 108 w1_2=w1_1;w1_1=w1; 109 b1_2=b1_1;b1_1=b1; 110 w2_2=w2_1;w2_1=w2; 111 b2_2=b2_1;b2_1=b2; 112 end 113 end 114 115 116 %% 产生双螺旋测试数据 117 %% 产生双螺旋数据,每类100个样本点，共200个样本 118 i=(1.5:1:test_number+0.5)'; %每类51个样本 119 120 %双螺旋数据点的产生方程 121 alpha2=pi*(i-1)/25; 122 beta2=0.4*((105-i)/104); 123 m0=0.5+beta2.*sin(alpha2); 124 n0=0.5+beta2.*cos(alpha2); 125 s0=zeros(test_number,1); 126 m1=0.5-beta2.*sin(alpha2); 127 n1=0.5-beta2.*cos(alpha2); 128 s1=ones(test_number,1); 129 130 test=[m0 n0 s0;m1,n1,s1]; %1条螺旋线数据点,3*102的矩阵 131 test_label1=test(:,end)'; %测试数据类别，1*102的行向量 132 test_data1=test(:,1:end-1)'; %测试数据属性，2*102的矩阵 133 134 %把1维的输出变成2维的输出,train_labe2为200*2的矩阵 135 for i=1:2*test_number 136 switch test_label1(i) 137 case 0 138 test_label2(i,:)=[1 0]; 139 case 1 140 test_label2(i,:)=[0 1]; 141 end 142 end 143 144 test_label=test_label2'; %test_label为2*102的矩阵 145 146 %% 画出测试数据双螺旋曲线 147 plot(m0,n0,'c+'); 148 hold on; 149 plot(m1,n1,'yo'); 150 legend('训练数据螺旋线1','训练数据螺旋线2','测试数据螺旋线1','测试数据螺旋线2'); 151 152 test_data=mapminmax('apply',test_data1,train_datas); 153 154 % %% 用训练到的模型对训练数据本身进行预测 155 % for i=1:102 156 % for j=1:midnum 157 % I(j)=train_data(:,i)'*w1(j,:)'+b1(j); 158 % Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量 159 % end 160 % predict(:,i)=w2'*Iout'+b2;%predict为2*102的矩阵 161 % end 162 % 163 % test_data=mapminmax('apply',train_data1,train_datas); 164 % test_label=train_label; 165 % test_label1=trian_label1; 166 167 %% 用训练到的模型预测数据 168 for i=1:2*test_number 169 for j=1:midnum 170 I(j)=test_data(:,i)'*w1(j,:)'+b1(j); 171 Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量 172 end 173 predict(:,i)=w2'*Iout'+b2;%predict为2*102的矩阵 174 end 175 176 %% 预测结果分析 177 for i=1:2*test_number 178 output_pred(i)=find(predict(:,i)==max(predict(:,i))); %out_pred为1*102的矩阵 179 end 180 181 error=output_pred-test_label1-1; % 182 183 184 %% 计算出每一类预测错误的个数总和 185 k=zeros(1,2); %k=[0 0] 186 for i=1:2*test_number 187 if error(i)~=0 %matlab中不能用if error(i)！=0 188 [b c]=max(test_label(:,i)); 189 switch c 190 case 1 191 k(1)=k(1)+1; 192 case 2 193 k(2)=k(2)+1; 194 end 195 end 196 end 197 198 199 %% 求出每一类总体的个数和 200 kk=zeros(1,2); %k=[0 0] 201 for i=1:2*test_number 202 [b c]=max(test_label(:,i)); 203 switch c 204 case 1 205 kk(1)=kk(1)+1; 206 case 2 207 kk(2)=kk(2)+1; 208 end 209 end 210 211 212 %% 计算每一类的正确率 213 accuracy=(kk-k)./kk 214
实验总结：
可以看出，自己实现的BP网络其预测效果教matlab自带的差些，这主要是其中一些算法没有进行优化，比如说初始权值的选取，自己实现时是随机产生的，这样的话不同的权值在相同迭代次数的情况下效果显然不同。另外本次试验的训练样本数才200个，也比较少，且中间层神经网络只有5个，进一步影响了试验的效果。
所以下一步的工作是，不管是matlab自带的工具还是自己写的函数，对参数的寻优过程要做一定的优化。
参考文献：
- W Zhao, DS Huang，The structure optimization of radial basis probabilistic neural networks based on genetic algorithms.
- Matlab中文论坛，MATLAB神经网络30个案例分析。
大家有更好的解法欢迎提出并交流。