rbf.m

matlab,神经网络

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

clc;
clear 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.5)*randn(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.5)*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 = 200;   % 最大学习次数
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 = 0.5*((Y-P_M1*W).^2);
E_sum1(1) = sum(E);

% 粗调阶段
k = 1;               % 隐节点数
k_Hide(k) = k;
while C_A <= C_max          % 判断粗调是否停止
    k = k+1;
    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 = 0.5*((Y-P_M1*W).^2);
    E_sum1(k) = sum(E);       % 当前训练样本的总误差
end
W = inv(P_M1'*P_M1+lambda*eye(size(P_M1'*P_M1)))*P_M1'*Y;   % 最优权值
E = 0.5*((Y-P_M1*W).^2);
E_sum(1) = sum(E);       % 当前训练样本的总误差
A1(1) = 0;                     % 当前的平均误差

figure;
plot(k_Hide,E_sum1);
xlabel('隐节点数'); ylabel('学习总误差'); title('粗调阶段学习总误差与隐节点数的关系');
xlim([0 k]);

%=============开始精调================
for Train_N = 2:Train_Nmax 
    clear P_M; clear P_M1; clear P_M11;
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*eye(size(P_M1'*P_M1)))*P_M1'*Y;    % 精调后的权值

% 计算当前所有训练样本的总误差和平均误差
E = 0.5*((Y-P_M1*W).^2);
E_sum(Train_N) = sum(E);       % 当前训练样本的总误差
A1(Train_N) = miu*A1(Train_N-1)+(1-miu)*E_sum(Train_N);    % 当前训练样本的平均误差

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

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

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

%====================测试网络=====================
for i = 1:Sample_N
    for j = 1:k
        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;
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('无噪声时的期望输出','有噪声时的期望输出','测试样本的实际输出');