一、RBF简介
1 RBF 1.1 什么是径向基函数 1985年，Powell提出了多变量插值的径向基函数（RBF）方法。径向基函数是一个取值仅仅依赖于离原点距离的实值函数，也就是Φ（x）=Φ(‖x‖),或者还可以是到任意一点c的距离，c点称为中心点，也就是Φ（x，c）=Φ(‖x-c‖)。任意一个满足Φ（x）=Φ(‖x‖)特性的函数Φ都叫做径向基函数，标准的一般使用欧氏距离（也叫做欧式径向基函数），尽管其他距离函数也是可以的。最常用的径向基函数是高斯核函数 ,形式为 k(||x-xc||)=exp{- ||x-xc||^2/(2*σ)^2) } 其中x_c为核函数中心,σ为函数的宽度参数 , 控制了函数的径向作用范围。
1.2 RBF神经网络
RBF神将网络是一种三层神经网络，其包括输入层、隐层、输出层。从输入空间到隐层空间的变换是非线性的，而从隐层空间到输出层空间变换是线性的。流图如下：
RBF网络的基本思想是：用RBF作为隐单元的“基”构成隐含层空间，这样就可以将输入矢量直接映射到隐空间，而不需要通过权连接。当RBF的中心点确定以后，这种映射关系也就确定了。而隐含层空间到输出空间的映射是线性的，即网络的输出是隐单元输出的线性加权和，此处的权即为网络可调参数。其中，隐含层的作用是把向量从低维度的p映射到高维度的h，这样低维度线性不可分的情况到高维度就可以变得线性可分了，主要就是核函数的思想。这样，网络由输入到输出的映射是非线性的，而网络输出对可调参数而言却又是线性的。网络的权就可由线性方程组直接解出，从而大大加快学习速度并避免局部极小问题。
径向基神经网络的激活函数可表示为：
其中xp为第p个输入样本，ci为第i个中心点，h为隐含层的结点数，n是输出的样本数或分类数。径向基神经网络的结构可得到网络的输出为：
当然，采用最小二乘的损失函数表示：
3 RBF神经网络的学习问题
求解的参数有3个：基函数的中心、方差以及隐含层到输出层的权值。 (1）自组织选取中心学习方法： 第一步：无监督学习过程，求解隐含层基函数的中心与方差 第二步：有监督学习过程，求解隐含层到输出层之间的权值 首先，选取h个中心做k-means聚类，对于高斯核函数的径向基，方差由公式求解： cmax为所选取中心点之间的最大距离。 隐含层至输出层之间的神经元的连接权值可以用最小二乘法直接计算得到，即对损失函数求解关于w的偏导数，使其等于0，可以化简得到计算公式为： （2）直接计算法 隐含层神经元的中心是随机地在输入样本中选取，且中心固定。一旦中心固定下来，隐含层神经元的输出便是已知的，这样的神经网络的连接权就可以通过求解线性方程组来确定。适用于样本数据的分布具有明显代表性。 （3）有监督学习算法 通过训练样本集来获得满足监督要求的网络中心和其他权重参数，经历一个误差修正学习的过程，与BP网络的学习原理一样，同样采用梯度下降法。因此RBF同样可以被当作BP神经网络的一种。
1.4 RBF神经网络与BP神经网络之间的区别 1.4.1 局部逼近与全局逼近 BP神经网络的隐节点采用输入模式与权向量的内积作为激活函数的自变量，而激活函数采用Sigmoid函数。各调参数对BP网络的输出具有同等地位的影响，因此BP神经网络是对非线性映射的全局逼近。 RBF神经网络的隐节点采用输入模式与中心向量的距离（如欧式距离）作为函数的自变量，并使用径向基函数（如Gaussian函数）作为激活函数。神经元的输入离径向基函数中心越远，神经元的激活程度就越低（高斯函数）。RBF网络的输出与部分调参数有关，譬如，一个wij值只影响一个yi的输出，RBF神经网络因此具有“局部映射”特性。
所谓局部逼近是指目标函数的逼近仅仅根据查询点附近的数据。而事实上，对于径向基网络，通常使用的是高斯径向基函数，函数图象是两边衰减且径向对称的，当选取的中心与查询点（即输入数据）很接近的时候才对输入有真正的映射作用，若中心与查询点很远的时候，欧式距离太大的情况下，输出的结果趋于0，所以真正起作用的点还是与查询点很近的点，所以是局部逼近；而BP网络对目标函数的逼近跟所有数据都相关，而不仅仅来自查询点附近的数据。
1.4.2 中间层数的区别 BP神经网络可以有多个隐含层，但是RBF只有一个隐含层。
1.4.3 训练速度的区别 使用RBF的训练速度快，一方面是因为隐含层较少，另一方面，局部逼近可以简化计算量。对于一个输入x，只有部分神经元会有响应，其他的都近似为0，对应的w就不用调参了。
1.4.4 RBF网络是连续函数的最佳逼近，而BP网络不是。
二、部分源代码
%% Mackey Glass Time Series Prediction using Spatio-Temporal Radial Basis Function (RBF) Neural Network
% Author: Shujaat Khan, shujaat123@gmail.com
clc;
clear all;
close all;
%% Loading Time Series Data
% I generated a series x(t) for t = 0,1, . . . ,3000, using mackey glass series equation with the following configurations: b = 0.1, a = 0.2, Tau = 20, and the initial conditions x(t - Tau) = 0.
load Dataset/Data.mat
% Training and Test datasets
time_steps=2; % prediction of #time_steps forward value (for this simple architechture time_steps<=3)
% Training
start_of_series_tr=100;
end_of_series_tr=2500;
% Test
start_of_series_ts=2500;
end_of_series_ts=3000;
P_train=Data(start_of_series_tr:end_of_series_tr-time_steps,2); % Input Data
f_train=Data(start_of_series_tr+time_steps:end_of_series_tr,2); % Label Data (desired output values)
indt=Data(start_of_series_tr+time_steps:end_of_series_tr,1);% Time index
SNR = 30; % signal to noise ratio
f_train=awgn(f_train,SNR); % Adding white Gaussian noise
%% Simulation parameters
% Defining architechture of the RBF-NN
[m n] = size(P_train);% Dimensions of input data [m]-length of signal, [n]-number of elements in each input
order=2; % Number of past values used for the prediction of future value
n1 = 10; % Number of hidden layer neurons
% Tuning parameters for training
epoch=10; % simulation rounds (number of times the same data pass through the NN for training)
eta=5e-2; % Gradient Descent step-size (learning rate)
runs=10; % Number of Monte Carlo simulations
Iti=[]; % Initial mean square error (MSE)
% Graphics/Plot parameters
fsize=13; % Fontsize
lw=2; % line width size
%% Training Phase
for run=1:runs % Monte Carlos simulations loop
% spread and centers of the Gaussian kernel
[temp, c, beeta] = kmeans(P_train,n1); % K-means clustering
% c=[c awgn(c,10)];
beeta=2*beeta; % Increasing spread of Gaussian kernel
% Initialization of weights and bias
w=randn(order,n1); % weight
b=randn(); % bias
for k=1:epoch % simulation rounds loop
I(k)=0; % reset MSE
U=zeros(1,order); % reset input vector
for i1=1:m % Iteration loop
% sliding window (updating input vector)
U(1:end-1)=U(2:end);
U(end)=P_train(i1); % current value of time-series
% Gaussian Kernel
for i2=1:n1
phi(:,i2)=exp((-(abs(U-c(i2,:))))/beeta(i2,:).^2);
end
% Calculate output of the RBF
y_train(i1)=sum(sum(w.*phi))+b;
e(i1)=f_train(i1)-y_train(i1); % instantaneous error in the prediction
% Gradient descent-based weight-update rule
w=w+eta*e(i1)*phi;
b=b+eta*e(i1);
% Mean square error
I(i1)=mse(e(1:i1)); % Objective Function
end
Itti(epoch,:)=I; % MSE for all iterations
end
Iti(run,:)=mean(Itti,1); % Mean MSE for all epochs
end
It=mean(Iti,1); % Mean MSE for all independent runs (Monte Carlo simulations)
%% Test Phase
P_test=Data(start_of_series_ts:end_of_series_ts-time_steps,2);
f_test=Data(start_of_series_ts+time_steps:end_of_series_ts,2);
indts=Data(start_of_series_ts+time_steps:end_of_series_ts,1);
[m n] = size(P_test);
for i1=1:m % Iteration loop
% sliding window (updating input vector)
U(1:end-1)=U(2:end);
U(end)=P_test(i1);
for i2=1:n1
phi(:,i2)=exp((-(abs(U-c(i2,:))))/beeta(i2,:).^2);
end
y_test(i1)=sum(sum(w.*phi))+b;
e_test(i1)=real(f_test(i1)-y_test(i1));
I(2400+i1)=mse(e_test(1:i1));
end
save ST_RBF.mat
% %% Results
% % Input and output signals (training phase)
% figure
% plot(indt,f_train,'k','linewidth',lw);
% hold on;
% plot(indt,y_train,'.:b','linewidth',lw);
% xlim([start_of_series_tr+time_steps end_of_series_tr]);
% h=legend('Actual Value (Training)','RBF Predicted (Training)','Location','Best');
% grid minor
% xlabel('Sample #','FontSize',fsize);
% ylabel('Magnitude','FontSize',fsize);
% set(h,'FontSize',12)
% set(gca,'FontSize',13)
% saveas(gcf,strcat('Time_SeriesTraining.png'),'png')
%
% % Input and output signals (test phase)
% figure
% plot(indts,f_test,'k','linewidth',lw);
% hold on;
% plot(indts,y_test,'.:b','linewidth',lw);
% xlim([start_of_series_ts+time_steps end_of_series_ts]);
% h=legend('Actual Value (Testing)','RBF Predicted (Testing)','Location','Best');
% grid minor
% xlabel('Sample #','FontSize',fsize);
% ylabel('Magnitude','FontSize',fsize);
% set(h,'FontSize',12)
% set(gca,'FontSize',13)
% saveas(gcf,strcat('Time_SeriesTesting.png'),'png')
%
% % Objective function (MSE) (training phase)
% figure
% plot(start_of_series_tr:end_of_series_tr-1,10*log10(I(1:end_of_series_tr-start_of_series_tr)),'+-b','linewidth',lw)
% h=legend('RBF (Training)','Location','North');
% grid minor
% xlabel('Sample #','FontSize',fsize);
% ylabel('MSE (dB)','FontSize',fsize);
% set(h,'FontSize',12)
% set(gca,'FontSize',13)
% saveas(gcf,strcat('Time_SeriesTrainingMSE.png'),'png')
%
% % Objective function (MSE) (test phase)
% figure
% plot(start_of_series_ts+time_steps:end_of_series_ts,10*log10(I(end_of_series_tr-start_of_series_tr+1:end)),'.:b','linewidth',lw+1)
% h=legend('RBF (Testing)','Location','South');
% grid minor
% xlabel('Sample #','FontSize',fsize);
% ylabel('MSE (dB)','FontSize',fsize);
% set(h,'FontSize',12)
% set(gca,'FontSize',13)
% saveas(gcf,strcat('Time_SeriesTestingMSE.png'),'png')
%
% % Mean square error
% 10*log10(((f_train'-y_train)*(f_train'-y_train)')/length(y_train))
% 10*log10(((f_test'-y_test)*(f_test'-y_test)')/length(y_test))
Results_graphs
三、运行结果
四、matlab版本及参考文献
1 matlab版本 2014a
2 参考文献 [1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016. [2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017. [3]周品.MATLAB 神经网络设计与应用[M].清华大学出版社，2013. [4]陈明.MATLAB神经网络原理与实例精解[M].清华大学出版社，2013. [5]方清城.MATLAB R2016a神经网络设计与应用28个案例分析[M].清华大学出版社，2018.
