神经网络学习笔记_1(BP网络分类双螺旋线)

这个也是本学期模式识别的一个小作业题。

实验目的：采用神经网络的方法对二维双螺旋样本点进行分类。

实验数据：由于本次的实验数据是双螺旋，需要用数学公式产生，其产生方法见参考文献[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)：

%% 产生双螺旋数据

studytra_num=100;  %训练样本数目
predict_test_number=100;    %测试样本数目

i=(1:1:studytra_num)';
%双螺旋数据点的产生方程
alpha1=pi*(i-1)/25;
beta=0.4*((105-i)/104);
x0=0.5+beta.*sin(alpha1);
y0=0.5+beta.*cos(alpha1);
z0=zeros(studytra_num,1);
x1=0.5-beta.*sin(alpha1);
y1=0.5-beta.*cos(alpha1);
z1=ones(studytra_num,1);

%% 事实证明BP神经网络在训练数据时与输入数据正负样本的顺序是有关系的
% 如果一开始的一半数据都是正(负)样本，后面的全是负(正)样本，则训练出来的
%效果不好，所以这里需要随机打乱
k=rand(1,2*studytra_num);
[m,n]=sort(k);

studytra=[x0 y0 z0;x1,y1,z1]; %1条螺旋线数据点,200*3的矩阵
trian_label1=studytra(n(1:2*studytra_num),end)';   %训练数据类别，1*200的行向量
studytra_data1=studytra(n(1:2*studytra_num),1:end-1)'; %训练数据属性，2*200的矩阵

%把1维的输出变成2维的输出,studytra_labe2为200*2的矩阵
for i=1:2*studytra_num
    switch trian_label1(i)
        case 0
            studytra_label2(i,:)=[1 0];
        case 1
            studytra_label2(i,:)=[0 1];
    end
end

[studytra_data,studytra_datas]=mapminmax(studytra_data1);%studytra_data为2*200的矩阵
studytra_label=studytra_label2'; %studytra_label为2*200的矩阵
         
plot(x0,y0,'r+');
hold on;
plot(x1,y1,'go');


%% 产生双螺旋测试数据

i=(1.5:1:predict_test_number+0.5)';    %每类51个样本
%双螺旋数据点的产生方程
alpha2=pi*(i-1)/25;
beta2=0.4*((105-i)/104);
m0=0.5+beta2.*sin(alpha2);
n0=0.5+beta2.*cos(alpha2);
s0=zeros(predict_test_number,1);
m1=0.5-beta2.*sin(alpha2);
n1=0.5-beta2.*cos(alpha2);
s1=ones(predict_test_number,1);

predict_test=[m0 n0 s0;m1,n1,s1]; %1条螺旋线数据点,3*102的矩阵
predict_test_label1=predict_test(:,end)';   %测试数据类别，1*102的行向量
predict_test_data1=predict_test(:,1:end-1)'; %测试数据属性，2*102的矩阵

%把1维的输出变成2维的输出,studytra_labe2为200*2的矩阵
for i=1:2*predict_test_number
    switch predict_test_label1(i)
        case 0
            predict_test_label2(i,:)=[1 0];
        case 1
            predict_test_label2(i,:)=[0 1];
    end
end

predict_test_label=predict_test_label2'; %predict_test_label为2*102的矩阵
         
%%  画出测试数据双螺旋曲线
plot(m0,n0,'c+');
hold on;
plot(m1,n1,'yo');
legend('训练数据螺旋线1','训练数据螺旋线2','测试数据螺旋线1','测试数据螺旋线2');

predict_test_data=mapminmax('apply',predict_test_data1,studytra_datas);

%% 至此，训练和测试数据如下所示：
studytra_data  ;    %2*200的矩阵
studytra_label  ;   %2*200的矩阵
predict_test_data   ;    %2*200的矩阵
predict_test_label   ;   %2*200的矩阵

% %%  调用matlab的神经网络工具箱函数用来训练,方法1
% net=newff(minmax(studytra_data),[10,2],{'tansig','purelin','studytragdm'})
% inputWeights=net.IW{1,1};
% inputbias=net.b{1};
% layerWeights=net.LW{2,1};
% layerbias=net.b{2};
% 
% net.studytraParam.show=50;
% net.studytraParam.lr=0.05;
% net.studytraParam.mc=0.9;
% net.studytraParam.epochs=1000;
% net.studytraParam.goal=1e-2;
% 
% [net,tr]=studytra(net,studytra_data,studytra_label);


%% 调用matlab的神经网络工具箱函数用来训练,方法2
net=feedforwardnet(8);
%注意此时进入train函数的样本每一列为一个样本的属性，列数为样本数，且命名时最好不要含有train,test字样，否则matlab会报错
net=train(net,studytra_data,studytra_label);
%view(net);
predict_label=sim(net,predict_test_data);


% %% 用训练到的模型预测数据
% for i=1:2*predict_test_number
%     for j=1:midnum
%         I(j)=predict_test_data(:,i)'*w1(j,:)'+b1(j);
%         Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量
%     end
%     predict_test(:,i)=w2'*Iout'+b2;%predict_test为2*102的矩阵
% end


%% 预测结果分析
for i=1:2*predict_test_number
    output_pred(i)=find(predict_label(:,i)==max(predict_label(:,i)));    %out_pred为1*102的矩阵
end

error=output_pred-predict_test_label1-1;    %


%% 计算出每一类预测错误的个数总和
k=zeros(1,2); %k=[0 0]
for i=1:2*predict_test_number
    if error(i)~=0    %matlab中不能用if error(i)！=0 
        [b c]=max(predict_test_label(:,i));
        switch c
            case 1
                k(1)=k(1)+1;
            case 2
                k(2)=k(2)+1;
        end
    end
end


%% 求出每一类总体的个数和
kk=zeros(1,2); %k=[0 0]
for i=1:2*predict_test_number
    [b c]=max(predict_test_label(:,i));
    switch c
        case 1
            kk(1)=kk(1)+1;
        case 2
            kk(2)=kk(2)+1;
    end
end


%% 计算每一类的预测正确率
accuracy=(kk-k)./kk

 

方法二：自己实现BP神经网络的训练和测试函数

 

实验说明：实现BP神经网络的方法参考的是文献[2]，主要分为下面6个步骤：

1. 网络的初始化
2. 隐含层输出计算
3. 输出层输出计算
4. 误差计算
5. 权值更新
6. 阈值更新

　　同样，训练和测试样本都为2*200的矩阵，且训练过程中设置的循环次数为200次，只有1个隐含层，且神经元个数为5个。

　　实验结果：

　　实验预测的精度如下所示：

　　

　　即第1条螺旋线的预测精度为94%,，第2条螺旋线的预测精度为37%。

实验代码(matlab)：

%% 产生双螺旋数据,每类100个样本点，共200个样本
train_num=100;
train_circle_number=5000;
test_number=100;
i=(1:1:train_num)';

%双螺旋数据点的产生方程
alpha1=pi*(i-1)/25;
beta=0.4*((105-i)/104);
x0=0.5+beta.*sin(alpha1);
y0=0.5+beta.*cos(alpha1);
z0=zeros(train_num,1);
x1=0.5-beta.*sin(alpha1);
y1=0.5-beta.*cos(alpha1);
z1=ones(train_num,1);

%% 事实证明BP神经网络在训练数据时与输入数据正负样本的顺序是有关系的
% 如果一开始的一半数据都是正(负)样本，后面的全是负(正)样本，则训练出来的
%效果不好，所以这里需要随机打乱
k=rand(1,2*train_num);
[m,n]=sort(k);

train=[x0 y0 z0;x1,y1,z1]; %1条螺旋线数据点,200*3的矩阵
trian_label1=train(n(1:2*train_num),end)';   %训练数据类别，1*200的行向量
train_data1=train(n(1:2*train_num),1:end-1)'; %训练数据属性，2*200的矩阵

%把1维的输出变成2维的输出,train_labe2为200*2的矩阵
for i=1:2*train_num
    switch trian_label1(i)
        case 0
            train_label2(i,:)=[1 0];
        case 1
            train_label2(i,:)=[0 1];
    end
end

train_label=train_label2'; %train_label为2*200的矩阵
         
plot(x0,y0,'r+');
hold on;
plot(x1,y1,'go');
%legend();

%% BP神经网络结构的初始化
%网络结构，2个输入，3个神经元，2个输出
innum=2;
midnum=5;
outnum=2;

[train_data,train_datas]=mapminmax(train_data1);

%输入输出取值阈值随机初始化
%w1矩阵表示每一行为一个隐含层神经元的输入权值
w1=rands(midnum,innum);%rands函数用来初始化神经元的权值和阈值是很合适的,w1为3*2的矩阵
b1=rands(midnum,1);%b1为3*1的矩阵
%w2矩阵表示每一列为一个输出层神经元的输入权值
w2=rands(midnum,outnum);%w2为3*2的矩阵
b2=rands(outnum,1);%b2为2*1的矩阵

%用来保存上一次的权值和阈值，因为后面的更新方差是递归的，要用到
w1_1=w1;w1_2=w1_1;
b1_1=b1;b1_2=b1_1;
w2_1=w2;w2_2=w2_1;
b2_1=b2;b2_2=b2_1;

%学习率的设定
alpha=0.05;

%训练10次就ok了，而不管训练后的结果如何
for train_circle=1:train_circle_number  ;
    for i=1:2*train_num; %200个训练样本
       %% 输入层的输出
        x=train_data(:,i);%取出第i个样本，x(i)为2*1的列向量
        %% 隐含层的输出
        for j=1:midnum;
            I(j)=train_data(:,i)'*w1(j,:)'+b1(j);  %I(j)为1*1的实数
            Iout(j)=1/(1+exp(-I(j)));   %Iout(j)也为1*1的实数
        end     %Iout为1*3的行向量   
        %% 输出层的输出
         yn=(Iout*w2)'+b2;   %yn为2*1的列向量，因此此时的传函为线性的，所以可以一步到位，不必上面
        
        %% 计算误差
        e=train_label(:,i)-yn; %e为2*1的列向量，保存的是误差值
        
        %计算权值变换率
        dw2=e*Iout; %dw2为2*3的矩阵，每一行表示输出接点的输入权值变化率
        db2=e'; %e为1*2的行向量
        
        for j=1:midnum
            S=1/(1+exp(-I(j)));
            FI(j)=S*(1-S);  %FI(j)为一实数，FI为1*3的行向量
        end
        
        for k=1:1:innum
            for j=1:midnum
                dw1(k,j)=FI(j)*x(k)*(e(1)*w2(j,1)+e(2)*w2(j,2));    %dw1为2*3的矩阵
                db1(j)=FI(j)*(e(1)*w2(j,1)+e(2)*w2(j,2));   %db1为1*3的矩阵
            end
        end
        
        %% 权值更新方程
        w1=w1_1+alpha*dw1'; %w1仍为3*2的矩阵
        b1=b1_1+alpha*db1'; %b1仍为3*1的矩阵
        w2=w2_1+alpha*dw2'; %w2仍为3*2的矩阵
        b2=b2_1+alpha*db2'; %b2仍为2*1的矩阵
        
        %% 保存上一次的权值和阈值
        w1_2=w1_1;w1_1=w1;
        b1_2=b1_1;b1_1=b1;
        w2_2=w2_1;w2_1=w2;
        b2_2=b2_1;b2_1=b2;
    end
end


%% 产生双螺旋测试数据
%% 产生双螺旋数据,每类100个样本点，共200个样本
i=(1.5:1:test_number+0.5)';    %每类51个样本

%双螺旋数据点的产生方程
alpha2=pi*(i-1)/25;
beta2=0.4*((105-i)/104);
m0=0.5+beta2.*sin(alpha2);
n0=0.5+beta2.*cos(alpha2);
s0=zeros(test_number,1);
m1=0.5-beta2.*sin(alpha2);
n1=0.5-beta2.*cos(alpha2);
s1=ones(test_number,1);

test=[m0 n0 s0;m1,n1,s1]; %1条螺旋线数据点,3*102的矩阵
test_label1=test(:,end)';   %测试数据类别，1*102的行向量
test_data1=test(:,1:end-1)'; %测试数据属性，2*102的矩阵

%把1维的输出变成2维的输出,train_labe2为200*2的矩阵
for i=1:2*test_number
    switch test_label1(i)
        case 0
            test_label2(i,:)=[1 0];
        case 1
            test_label2(i,:)=[0 1];
    end
end

test_label=test_label2'; %test_label为2*102的矩阵
         
%%  画出测试数据双螺旋曲线
plot(m0,n0,'c+');
hold on;
plot(m1,n1,'yo');
legend('训练数据螺旋线1','训练数据螺旋线2','测试数据螺旋线1','测试数据螺旋线2');

test_data=mapminmax('apply',test_data1,train_datas);

% %% 用训练到的模型对训练数据本身进行预测
% for i=1:102
%     for j=1:midnum
%         I(j)=train_data(:,i)'*w1(j,:)'+b1(j);
%         Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量
%     end
%     predict(:,i)=w2'*Iout'+b2;%predict为2*102的矩阵
% end
% 
% test_data=mapminmax('apply',train_data1,train_datas);
% test_label=train_label;
% test_label1=trian_label1;

%% 用训练到的模型预测数据
for i=1:2*test_number
    for j=1:midnum
        I(j)=test_data(:,i)'*w1(j,:)'+b1(j);
        Iout(j)=1/(1+exp(-I(j)));%Iout为1*3的行向量
    end
    predict(:,i)=w2'*Iout'+b2;%predict为2*102的矩阵
end

%% 预测结果分析
for i=1:2*test_number
    output_pred(i)=find(predict(:,i)==max(predict(:,i)));    %out_pred为1*102的矩阵
end

error=output_pred-test_label1-1;    %


%% 计算出每一类预测错误的个数总和
k=zeros(1,2); %k=[0 0]
for i=1:2*test_number
    if error(i)~=0    %matlab中不能用if error(i)！=0 
        [b c]=max(test_label(:,i));
        switch c
            case 1
                k(1)=k(1)+1;
            case 2
                k(2)=k(2)+1;
        end
    end
end


%% 求出每一类总体的个数和
kk=zeros(1,2); %k=[0 0]
for i=1:2*test_number
    [b c]=max(test_label(:,i));
    switch c
        case 1
            kk(1)=kk(1)+1;
        case 2
            kk(2)=kk(2)+1;
    end
end


%% 计算每一类的正确率
accuracy=(kk-k)./kk
               

 

实验总结：

　　可以看出，自己实现的BP网络其预测效果教matlab自带的差些，这主要是其中一些算法没有进行优化，比如说初始权值的选取，自己实现时是随机产生的，这样的话不同的权值在相同迭代次数的情况下效果显然不同。另外本次试验的训练样本数才200个，也比较少，且中间层神经网络只有5个，进一步影响了试验的效果。

　　所以下一步的工作是，不管是matlab自带的工具还是自己写的函数，对参数的寻优过程要做一定的优化。

 

参考文献：

1. W Zhao, DS Huang，The structure optimization of radial basis probabilistic neural networks based on genetic algorithms.
2. Matlab中文论坛，MATLAB神经网络30个案例分析。

大家有更好的解法欢迎提出并交流。