dynrbf.m

matlab,神经网络

%=====================RBF网络动态设计方法=====================

clc;
clear all;
close all;
% 产生训练样本
Sample_N = 100;    % 样本数
X = -4+8*rand(Sample_N,1);   % 服从区间(-4,4)内均匀分布的训练样本
Y1 = 1.1*(1-X+2*(X.^2)).*exp(-(X.^2)/2);  % 教师输出（无噪声）
Y = Y1+sqrt(0.25)*randn(Sample_N,1);      % 教师输出（有噪声）
% X1 = -4+8*rand(Sample_N,1);           % 测试样本
X1 = linspace(-4,4,100)';
Y_Ex = 1.1*(1-X1+2*(X1.^2)).*exp(-(X1.^2)/2);   % 期望输出（无噪声）
Y_Ex1 = Y_Ex+sqrt(0.25)*randn(Sample_N,1);     % 期望输出（有噪声）

% 设置参数
r = 1.5;   % 扩展常数
lambda = 3e-4;      % 正则化系数
C_max = 1e6;        % 条数极限
eta = 1e-4;         % 学习率
kk = 1.0;           % 重叠系数
w_min = 0.1;        % 临界权值
Train_Nmax = 150;   % 最大学习次数
D = 0;              % 期望误差
delta_lambda = 2e-3;   % 正则化系数增量
miu = 0.95;
rou = 0.95;

%===========开始训练=============
% 选择RBF网第一个数据中心，并计算输出权值
for j = 1:Sample_N
    for i = 1:Sample_N
        S(i,j) = exp(-norm(X(i)-X(j))/r);    % 设计阵
    end
end
E1 = Y'*S;      % Gaussian正则化方法时的目标函数
[a1 b1] = max(E1);
c(1) = X(b1);      % 选择第b（使得E最大）个输入样本x为RBF网络第一个数据中心
P_M(:,1) = S(:,b1);
P_M11 = P_M;
P_M11(:,end+1) = ones(Sample_N,1);
P_M11(:,end+1) = S(:,b1);
A = P_M11'*P_M11+lambda*eye(size(P_M11'*P_M11));
C_A = norm(A,'fro')*norm(inv(A),'fro');       % 条件数
P_M1 = P_M11(:,1:end-1);
W = inv(P_M1'*P_M1+lambda*eye(size(P_M1'*P_M1)))*P_M1'*Y;   % 最优权值
E(1) = (Y-P_M1*W)'*(Y-P_M1*W)+lambda*W(end-1)'*W(end-1);
A1(1) = (1-miu)*E(1);    % 当前训练样本的平均误差

%=====================开始粗调=======================
k = 1;               % 隐节点数
abc(k) = k;          % 训练次数
k_Hide(k) = k;
while C_A <= C_max          % 判断粗调是否停止
    k = k+1;
    abc(k) = k;
    k_Hide(k) = k;
    P_M1 = P_M;
    P_M1(:,end+1) = ones(Sample_N,1);
    for j = 1:Sample_N
        P_M11 = P_M1;
        P_M11(:,end+1) = S(:,j);
        H_M11 = eye(Sample_N)-P_M11*inv(P_M11'*P_M11+lambda*eye(size(P_M11'*P_M11)))*P_M11';
        E_M11(j) = Y'*H_M11*Y;
    end
    [a2 b2] = min(E_M11);
    c(k) = X(b2);      % 第k个数据中心    k = 2,3,4,...
    P_M(:,end+1) = S(:,b2);
    P_M11 = P_M;
    P_M11(:,end+1) = ones(Sample_N,1);
    P_M11(:,end+1) = S(:,b2);
%     A = P_M11'*P_M11;    
    A = P_M11'*P_M11+lambda*eye(size(P_M11'*P_M11));
    C_A = norm(A,'fro')*norm(inv(A),'fro');      % 条件数
    P_M1 = P_M11(:,1:end-1);
    W = inv(P_M1'*P_M1+lambda*eye(size(P_M1'*P_M1)))*P_M1'*Y;   % 权值
    E(abc(end)) = (Y-P_M1*W)'*(Y-P_M1*W)+lambda*W(end-1)'*W(end-1);      % 当前训练样本的总误差
    A1(abc(end)) = miu*A1(abc(end)-1)+(1-miu)*E(abc(end));    % 当前训练样本的平均误差
end
lambda(1:k) = lambda;
%=====================粗调结束=======================

%=====================开始精调=======================
for Train_N = 2:Train_Nmax-k+1
    clear P_M; clear P_M1; clear P_M11;
    abc(k+Train_N-1) = k+Train_N-1;
for i = 1:k
    % 调整数据中心
    k1 = 0;
    for j = 1:Sample_N
        if norm(X(j)-c(i),'fro')<kk*r
            k1 = k1+1;
            fx = W(i)*exp(-norm(X(j)-c(i),'fro')/r)+W(end);
            delta_c(k1) = 4*(eta/r)*(X(j)-c(i))*exp(-abs(Y(j)-fx-c(i))/r)*W(i);    % 数据中心调节量
        end
    end
    c(i) = c(i)+sum(delta_c);      % 精调之后的数据中心
    
    % 调整网络输出权值和输出偏移
    for j = 1:Sample_N
        P_M(j,i) = exp(-norm(X(j)-c(i),'fro')/r);
    end
    P_M1 = P_M;
    P_M1(:,end+1) = ones(Sample_N,1);
end
clear W;
W = inv(P_M1'*P_M1+lambda(abc(end)-1)*eye(size(P_M1'*P_M1)))*P_M1'*Y;    % 精调后的权值

% 对隐节点进行剪枝
k10 = 0;
for i = 1:k
    if abs(W(i))<w_min
        W(i) = 0;
        k10 = k10+1;
    end
end 
k = k-k10;         % 剪枝后的隐接点个数
k_Hide(abc(end)) = k;
clear c1; clear W1; clear P_M2;
c1 = c;
W1 = W;
P_M2 = P_M1;
o = 0;
for j = 1:length(c)               % 提取权值和数据中心中的非零值
    if W(j) == 0
        W1(j-o,:) = [];
        c1(j-o) = []; 
        P_M2(:,j-o) = [];
        o = o+1;
    end
end
clear c; clear W; clear P_M1;
c = c1;
W = W1;
P_M1 = P_M2;

% 计算当前所有训练样本的总误差和平均误差
E(abc(end)) = (Y-P_M1*W)'*(Y-P_M1*W)+lambda(abc(end)-1)*W(end-1)'*W(end-1);          % 当前训练样本的总误差
A1(abc(end)) = miu*A1(abc(end)-1)+(1-miu)*E(abc(end));    % 当前训练样本的平均误差

% 判断算法是否结束
if E(abc(end))==D
    break;
end

% 调整正则化系数
if E(abc(end))<E(abc(end)-1) || E(abc(end))<D
    lambda(abc(end)) = lambda(abc(end)-1)+delta_lambda;
elseif E(abc(end))>=E(abc(end)-1) && E(abc(end))<A1(abc(end)) && E(abc(end))>=D
    lambda(abc(end)) = lambda(abc(end)-1)-delta_lambda;
elseif E(abc(end))>=E(abc(end)-1) && E(abc(end))>=A1(abc(end)) && E(abc(end))>=D
    lambda(abc(end)) = rou*lambda(abc(end)-1);
end
end
%====================精调结束=====================

%====================测试网络=====================
for j = 1:k
    for i = 1:Sample_N
        P_M(i,j) = exp(-norm(X1(i)-c(j))/r);
    end
end
P_M1 = P_M;
P_M1(:,end+1) = ones(Sample_N,1);
Y_Test = P_M1*W;     % 测试样本的实际输出

%==================画图====================
figure
[s t] = max(k_Hide);
plot(abc,k_Hide);
xlabel('训练次数（N）'); ylabel('隐节点个数');
title('隐节点个数与训练次数的关系');
ylim([0 s+5]);
plotxline(abc(t),'r');grid on;

figure
plot(abc,lambda);
xlabel('训练次数（N）'); ylabel('正则化系数');
title('正则化系数与训练次数的关系');
ylim([0 0.08]);
plotxline(abc(t),'r');grid on;

% figure;
% plot(abc,E);
% xlabel('学习次数'); ylabel('学习误差'); title('学习误差与学习次数的关系');
% xlim([0 abc(end)]);
% plotxline(abc(t),'r');grid on;

figure
plot(abc,E,'b',abc,A1,'r');
xlabel('学习次数'); ylabel('误差'); title('误差与学习次数的关系');
xlim([0 abc(end)]);
plotxline(abc(t),'r');grid on;
legend('学习误差','平均误差');

figure;
hold on;grid on;
plot(X1,Y_Ex,'b');
plot(X1,Y_Ex1,'y');
plot(X1,Y_Test,'r');
xlabel('测试样本（X）'); ylabel('函数输出（F(X)）');
legend('无噪声时的期望输出','有噪声时的期望输出','测试样本的实际输出');title('DYNRBF拟合曲线');

disp('测试样本误差（期望输出无噪声时）---------->:');
disp((Y_Ex-Y_Test)'*(Y_Ex-Y_Test)+lambda(end)*W(end-1)'*W(end-1));
disp('测试样本误差（期望输出有噪声时）---------->:');
disp((Y_Ex1-Y_Test)'*(Y_Ex1-Y_Test)+lambda(end)*W(end-1)'*W(end-1));
disp('训练结束时的训练误差---------->:');
disp(E(end));