BP的matlab实现

%2015.04.26 Kang Yongxin ----v 2.0
%完成作业中BP算法，采用批量方式更新权重
%%
%输入数据格式
%x 矩阵 : 样本个数*特征维度 
%y 矩阵 ：样本个数*类别个数（用01000形式表示）
close all;
clear all;
clc ;
load data.mat;%
x_test=x(1:3:30,:);%从原始数据中留出一部分 作为测试样本
y_test=y(1:3:30,:);
x_train=[x(2:3:30,:);x(3:3:30,:)];%x(1:2:30,:);%[x(11:25,:);x(26:30,:);x(1:10,:)];
y_train=[y(2:3:30,:);y(3:3:30,:)];%%[y(11:25,:);y(26:30,:);y(1:10,:)];
%%
%定义变量名称,初始化网络
d=size(x_train,2);%特征维度，也是输入层节点个数
num_trains=size(x_train,1);%训练样本个数
n_class=size(y_train,2);%样本类别数
node_layer=[d 4 n_class];%每层的节点个数 %[d 3 5 n_class];%构建更多层的网络
num_layer=size(node_layer,2)-1;%网络层数
for i=1:1:num_layer-1
    f_name{i}='sigmoid';%对应每层的激活函数
end
f_name{num_layer}='tanh';%'sigmoid';%最后一层的激活函数
eta=0.08;%学习率
theta=10e-4;%终止条件
W=cell(num_layer,1);%初始化权重矩阵，都设曾1,
for i=1:1:num_layer
    W{i}=rand(node_layer(i),node_layer(i+1));
end
W_init=W;
%开始循环
item=1;
while item>0 && item<1500
  %%
  %初始化权值增量
   for layer=1:1:num_layer
        delta_sum{layer}=zeros(size(W{layer})) ;
   end
   %%
    for k=1:1:num_trains
        %对于每个样本的循环
       %%%%%%%%%%%%%%%%%%%
        %前向计算
        x_in=x_train(k,:);
        for layer=1:1:num_layer
            %对于每一层进行前向计算，并且保存输出值y_out
            y_out{layer}=forward(x_in,W{layer},f_name{layer});%批量更新时候W要跟随外层循环
            x_in=y_out{layer};
        end
        %%%%%%%%%%%%%%%
        %%%%%%%%%%%%%% 
        %反向传播，最后一层要单算，因为只有一个输出 
        delta_out{layer}=y_train(k,:)-y_out{layer};%输出与真值的差
        J(k)=0.5*sum(delta_out{layer}.^2);%均方误
        delta_error{layer}=delta_out{layer}.*d_function(y_out{layer},f_name{num_layer});%从指向节点收集到的误差
        delta_w{layer}=eta*(y_out{layer-1})'*delta_error{layer};%本层的权重变化量
        while layer>1
           %反向传播误差，保存delta_w
           layer=layer-1;
           delta_error{layer}=delta_error{layer+1}*(W{layer+1})'.*d_function(y_out{layer},f_name{layer});%从指向节点收集到的误差并使用之前的权重进行加权
           if layer~=1
               %如果没到第一层，就用layer-1层的输出作为输入
               delta_w{layer}=eta*(y_out{layer-1})'*delta_error{layer};%本层的权重变化量
           else
               %如果是第一层就用本次训练的x作为输入   
               delta_w{layer}=eta*(x_train(k,:))'*delta_error{layer};%本层的权重变化量
           end
        end  
       %%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%
       %批量更新，要对所有样的贡献进行加和
       for layer=1:1:num_layer
            delta_sum{layer}=delta_sum{layer}+delta_w{layer} ;
       end
    end%k个样本对权重 变化量的贡献计算完毕
    %%
    
    figure(1);
    JW(item)=sum(J);
    if item>10
        Delta_JW=abs(JW(item)-JW(item-1));
        if Delta_JW<theta
             break;%循环终止条件
        end
    end
    plot(item,JW(item),'.')
    hold on;
    item=item+1;
    
  %更新权重
    for layer=1:1:num_layer
        %对于每一层更新权重,这个更新放的位置决定是进行批量更新还是每次更新
       W{layer}=W{layer}+delta_sum{layer};%批量更新
    end
end
%%
%计算准确率
x_in=x_test;
for layer=1:1:num_layer
    %对于每一层进行前向计算，并且保存输出值y_out
    y_out{layer}=forward(x_in,W{layer},f_name{layer});%批量更新时候W要跟随外层循环
    x_in=y_out{layer};
end
[C,I]=max(y_out{layer},[],2);        
[C,I_true]=max(y_test,[],2);
Trues=find((I-I_true)==0);
Precision=size(Trues,1)/size(y_test,1)

function [ y ] = forward( x,w,f_name )
%FORWARD 前向计算得到输出值
%   输入 x:输入向量1*m
%        w:权重矩阵m*n
%        f_name:函数名称（暂时支持'sigmoid''tanh'）
%   输出 y:输出向量1*n 公式为  y=f(x*w)
if strcmp(f_name,'sigmoid')
    y=sigmoid(x*w);
else if strcmp(f_name,'tanh')
        y=tanh(x*w);%matlab 自带
    else
       disp('wrong function name ');
    end
end

end

function [d_f ] = d_function( y,f_name )
%D_FOUNCTION 对相关函数进行求导数
%   输入是该层的输出值y 
if strcmp(f_name,'sigmoid')
    d_f=y.*(1-y);
else if strcmp(f_name,'tanh')
        d_f=1-y.^2;%matlab 自带
    else
       disp('wrong function name at d_function');
    end
end
end

function [ y ] = sigmoid( x )
%MY_SIGMOID 
%   输入向量x，输出
y=1./(1+exp(-x));
end