BP及RBF神经网络代码及时间序列预测问题

神经网络第一次作业130337杨康题目二：MG时间序列预测问题一、题目描述用MLP或者RBF网实现函数逼近，Mackey-Glass时间序列预测建模。时间序列由以下时延微分方程产生：dx(t)dt=0.2𝑥(𝑡−𝜏)1+𝑥10(𝑡−𝜏)−0.1𝑥(𝑡)令时延参数𝜏=17，初始条件𝑥(0)=1.2，且当𝑡0时x(t)=0。可以用4阶龙格库塔方法离散化。设计MLP或PBF网，使之能根据输入[x(t-18)x(t-12)x(t-6)x(t)]预测x(t+85)的值。用3000个样本进行训练（t=201到3200），用500个样本进行测试（t=5001到5500）。二、研究进展所谓时间序列，又称动态数据，指一组按时间顺序排列的数字序列，数据带有随机性，相互之间存在某种统计上的联系。时间序列按照性态可分为有确定规律的、混沌的和完全随机的。20多年来，国内外许多学者对时间序列的建模预测做了很多工作。90年代初以前，在数学界和工程界许多学者的共同努力下。国外以1976年GeorgeE.P.Box和GwilymM.Jenkins等的专著《TimeSeriesAnalysis:ForecastingandControl》和1983年S.M.Pandit和ShienMingWu的专著《TimeSeriesandSystemAnalysiswithApplications》为标志。国内以1983年安鸿志、陈兆国的专著《时间序列的分析与应用》和1991年杨叔子等著的《时间序列分析的工程应用》等为标志。此外，还有专门的期刊JournalofTimeSeriesAnalysis、国际会议和专题讨论，使时间序列建模预测技术从理论到应用都已经有了长足的发展。Padgett提出神经网络在时间序列预测中几种典型应用；EsparciaAlcazar提出利用遗传算法优化多层神经网络结构，并将此神经网络用于时间序列预测；A.S.Pandya对有噪声的时间序列预测进行了研究；FrancescoMasulli研究了模糊神经系统，并利用该系统对时间序列进行预测；RomanRosipal利用资源分配型RBF神经网络对混沌时间序列进行预测；吴春国等人基于正规正交分解(ProperOrthogonalDecomposition，POD)提出一种适用于非线性时间序列预测的径向基函数（RBF)神经网络模型-POD-RBF神经网络模型。该模型在选取中心时考虑了时间序列数据之间的时序关系，并且使得中心的选取具有并行性。北京航空航天大学的柳萍等人针对混沌时间序列预测问题，提出了一种基子小波框架的小波核函数与最小二乘支持向量机相结合的方法。该方法不仅能够以较高的精度逼近任意函数，而且还适用于混沌信号的局部分析，提高了最小二乘支持向量机的模型泛化能力。南京航空航天大学的张军峰基于一种新型聚类算法的RBF神经网络对混沌时间序列进行预测。文中提出了一种基于高斯基的距离度量,并联合输入输出聚类的策略。基于Fisher可分离率设计高斯基距离度量中的惩罚因子，可以提高聚类的性能。而输入输出聚类策略的引入，建立了聚类性能与网络预测性能之间的联系。三、算法及代码3.1竞争BP算法算法流程：（1）初始化权值；（2）依次输入3000个学习样本；（3）依次计算出各层的输出；（4）求出各层的反传误差；（5）比较隐层各单元的误差，最大的误差对应的权矢量正常修正，其他神经元的权矢量都与最大单元相反的方向修正；（6）按权值修正公式修正各层的权值和阈值；（7）按照新的权值在此计算，直至误差满足要求或者达到最大训练次数；（8）将测试的500个样本输入进去，得到测试结果，求出测试误差。程序代码：functionMackeyGlassN=60000;t=zeros(N,1);x=zeros(N,1);a=0.2;b=0.1;x(1)=1.2;t(1)=0;h=0.1;tau=17;fork=1:N-1;t(k+1)=t(k)+h;ift(k)tau;k1=-b*x(k);k2=-b*(x(k)+h*k1/2);k3=-b*(x(k)+k2*h/2);k4=-b*(x(k)+k3*h);x(k+1)=x(k)+(k1+2*k2+2*k3+k4)*h/6;elsen=floor((t(k)-tau-t(1))/h+1);k1=f(x(n))-b*x(k);k2=f(x(n))-b*(x(k)+h*k1/2);k3=f(x(n))-b*(x(k)+k2*h/2);k4=f(x(n))-b*(x(k)+k3*h);x(k+1)=x(k)+(k1+2*k2+2*k3+k4)*h/6;endendfigureholdongridplot(t,x);SamNum=3000;h=0.1;TestSamNum=500;HiddenUnitNum=20;InDim=4;OutDim=1;rand('state',sum(3000*clock))NoiseVar=0.1;Noise=NoiseVar*rand(OutDim,SamNum);SamIn=rand(InDim,SamNum);SamOutNoNoise=rand(OutDim,SamNum);TestSamIn=rand(InDim,TestSamNum);TestSamOut=rand(OutDim,TestSamNum);fork=201:3200SamIn(1,k-200)=x((k-18)/h+1);SamIn(2,k-200)=x((k-12)/h+1);SamIn(3,k-200)=x((k-6)/h+1);SamIn(4,k-200)=x((k)/h+1);SamOutNoNoise(k-200)=x((k+85)/h+1);endSamOut=SamOutNoNoise+Noise;forl=5001:5500TestSamIn(1,l-5000)=x((l-18)/h+1);TestSamIn(2,l-5000)=x((l-12)/h+1);TestSamIn(3,l-5000)=x((l-6)/h+1);TestSamIn(4,l-5000)=x((l)/h+1);TestSamOut(l-5000)=x((l+85)/h+1);endMaxEpochs=20000;Ir=0.00005;E0=0.5;W1=0.1*rand(HiddenUnitNum,InDim);B1=0.1*rand(HiddenUnitNum,1);W2=0.1*rand(OutDim,HiddenUnitNum);B2=0.1*rand(OutDim,1);W1Ex=[W1B1];W2Ex=[W2B2];SamInEx=[SamIn'ones(SamNum,1)]';ErrHistory=[];fori=1:MaxEpochsHiddenOut=logsig(W1Ex*SamInEx);HiddenOutEx=[HiddenOut'ones(SamNum,1)]';NetworkOut=W2Ex*HiddenOutEx;Error=SamOut-NetworkOut;SSE=sumsqr(Error)ErrHistory=[ErrHistorySSE];ifSSEE0,break,endDelta2=Error;Delta1=W2'*Delta2.*HiddenOut.*(1-HiddenOut);max1=0;d=1;forj=1:HiddenUnitNumnorm(Delta1(j,:),2);ifmax1=norm(Delta1(j,:),2)max1=norm(Delta1(j,:),2);d=j;endendforj=1:HiddenUnitNumifj==dDelta1(j,:)=Delta1(d,:);elseDelta1(j,:)=-Delta1(d,:)/4;endenddW2Ex=Delta2*HiddenOutEx';dW1Ex=Delta1*SamInEx';W1Ex=W1Ex+Ir*dW1Ex;W2Ex=W2Ex+Ir*dW2Ex;W2=W2Ex(:,1:HiddenUnitNum);endW1=W1Ex(:,1:InDim)B1=W1Ex(:,InDim+1)W2B2=W2Ex(:,1+HiddenUnitNum);TestHiddenOut=logsig(W1*TestSamIn+repmat(B1,1,TestSamNum));TestNNOut=W2*TestHiddenOut+repmat(B2,1,TestSamNum);figureholdongridplot(5001:5500,TestNNOut,'k-')plot(5001:5500,TestSamOut,'k--')figureholdongrid[xx,Num]=size(ErrHistory);plot(1:Num,ErrHistory,'k-');figureholdongridplot(5001:5500,TestSamOut-TestNNOut,'k-');BPTSSE=sqrt((sumsqr(TestSamOut-TestNNOut)/500))functiony=f(x)a=0.2;y=a*x/(1+x^10);说明：在竞争BP算法中，隐层数为1，单层隐节点数为20；最大训练次数为20000，学习率为0.0005，误差要求为0.5；在隐层权值修正中，非最大误差神经元权矢量修正为−δ/4；采用RSME测试误差。3.2基于聚类的RBF网络算法流程：（1）初始化：从样本输入中选取前88个样本作为初始聚类中心；（2）计算所有样本输入与聚类中心的距离；（3）对输入样本进行分类；（4）重新计算各类的新的聚类中心；（5）若该次聚类中心与上次重合则进入下一步；否则返回到步骤2；（6）确定各隐节点的扩展常数；（7）将测试的500个样本输入进去，得到测试结果，求出测试误差。程序代码：functionMackeyGlassN=60000;t=zeros(N,1);x=zeros(N,1);a=0.2;b=0.1;x(1)=1.2;t(1)=0;h=0.1;tau=17;fork=1:N-1;t(k+1)=t(k)+h;ift(k)tau;k1=-b*x(k);k2=-b*(x(k)+h*k1/2);k3=-b*(x(k)+k2*h/2);k4=-b*(x(k)+k3*h);x(k+1)=x(k)+(k1+2*k2+2*k3+k4)*h/6;elsen=floor((t(k)-tau-t(1))/h+1);k1=f(x(n))-b*x(k);k2=f(x(n))-b*(x(k)+h*k1/2);k3=f(x(n))-b*(x(k)+k2*h/2);k4=f(x(n))-b*(x(k)+k3*h);x(k+1)=x(k)+(k1+2*k2+2*k3+k4)*h/6;endendfigureholdongridplot(t,x);SamNum=3000;h=0.1;TestSamNum=500;InDim=4;ClusterNum=88;Overlap=1.0;rand('state',sum(3000*clock))NoiseVar=0.1;Noise=NoiseVar*rand(1,SamNum);SamIn=rand(InDim,SamNum);SamOutNoNoise=rand(1,SamNum);TestSamIn=rand(InDim,TestSamNum);TestSamOut=rand(1,TestSamNum);fork=201:3200SamIn(1,k-200)=x((k-18)/h+1);SamIn(2,k-200)=x((k-12)/h+1);SamIn(3,k-200)=x((k-6)/h+1);SamIn(4,k-200)=x((k)/h+1);SamOutNoNoise(k-200)=x((k+85)/h+1);endSam