hzh_nn_ar_emd_compose.m

利用神经网络预测时间序列

close all
%---------- Generate training and test set ----------
clear all
clc
load sunspot2003;%加载太阳黑子数
org_data=scale(sunspot2003,0,1); % 归一化load the normalised originate data sunspots
%-------------
load emd_sunspot
%x=imf;
[ro,vo]=size(imf); %得到域波分量imf的行数和列数
%%% -------------------*******------------------
for i=1:ro
  [x(i,:),xmin,xmax]=scale(imf(i,:),0,1);%归一化各个基本模式分量
end
        
title('Sunspot benchmark data after EMD decomposition');
for i=1:ro-1              %display the basic mode functions, leave out the residual component
  subplot(round((ro-1)/2),2,i);%共有 ro/2 行，每行显示两个分量
  plot(1:length(x),x(i,:));
  ylabel(['norm imf',num2str(i)]);
end

%=====================------------------------------
nlag=9;   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&时间序列的阶数 
train=cell(1,ro);  %Create cell array
test= cell(1,ro);
for i=1:ro   %对每个分量产生训练样本和检验样本
   [x0,trai,tes]=tsgenf_hzh(x(i,:),304,244,1,nlag);
   train{i}=trai;%每个分量的训练样本
   test{i}=tes;  %每个分量的检验样本
end
%------------******-----------------
%  W1    : Input-to-hidden layer weights. The matrix dimension is
%          dim(W1) = [(# of hidden units) * (inputs + 1)] (the 1 is due to the bias)
%  W2    : hidden-to-output layer weights.
%           dim(W2) = [(outputs)  *  (# of hidden units + 1)]
  rand('seed',0);
  
  NetDef = ['HHHHHH'       
            'L-----'];     
        
    maxiter = 1000;
    stop_crit = 1e-12;
    lambda=1;
    D=1;%加动量项有助于快速收敛
    trparms=[maxiter stop_crit lambda D];
           
     W1=cell(1,ro);
     W2=cell(1,ro);
   for i=1:ro 
     W1_in = rand(6,nlag+1)-0.5;  % Weights to hidden layer 权值加减-0。5 +0。5有助于预测的准确性
     W2_in = rand(1,7)-0.5;  % Weights to output
     [W1_out,W2_out,PI_vector,iter,lambda]=marq(NetDef,W1_in,W2_in,train{i}(:,1:nlag)',train{i}(:,nlag+1)',trparms);
     W1{i}=W1_out;
     W2{i}=W2_out;
     
     figure
    semilogy(PI_vector)
    title('Criterion evaluated after each iteration')
    xlabel('Iteration (epoch)')
    ylabel('Criterion')
    grid
    % -----------  Validate Network  -----------
    [Y_sim(i,:),E,PI] = nneval(NetDef,W1{i},W2{i},train{i}(:,1:nlag)',train{i}(:,nlag+1)',1);
    [Y_sim1(i,:),E,PI,tr,E1,PI1] = nneval_hzh(NetDef,W1{i},W2{i},test{i}(:,1:nlag)',test{i}(:,nlag+1)',1);
        
  end
  %  sum(E.*E)/length(train)
%    sum(E1.*E1)/length(test)


%-------------预测完所有分量后，对原始数据进行线性预测
net=newlind(Y_sim,org_data(nlag+1:244)');
 net_out1=sim(net,Y_sim);
figure;plot(net_out1,'r'),hold on
plot(org_data(nlag+1:244),'b'),hold off
e0=net_out1-org_data(nlag+1:244)';    
sum(e0.*e0)/length(net_out1)      %求训练样本的均方误差
%----------------------------------------------------------
pause
net_out2=sim(net,Y_sim1);
figure;plot(net_out2,'r'),hold on
plot(org_data(245:304),'b'),hold off

e=net_out2-org_data(245:304)';    %在244点的时候，性能超过单纯的神经网络，但比线性模型稍逊
sum(e.*e)/60                      %求检验样本的均方误差

%-------------------------------------------------------------------
net_out3=scale(net_out2,0,1);     %线性网络的输入归一化和结果显示
figure;plot(net_out3,'r'),hold on
plot(org_data(245:304),'b'),hold off
e1=net_out3-org_data(245:304)';
sum(e1.*e1)/60