NARX OPTIMUM HIDDEN NODES NUMBER

I used the code below to try to get the optimum number of hidden nodes for a narx network....I would like any suggestions concerning if I am implementing the code correctly. THANKS

        close all,clear all, clc, plt=0;

          tic

         load 'ProjectD'

        X = GasProduced; %[1x1500] cell
        T = OilRate;     %[1x1500] cell

         x  = cell2mat(X);
         t  = cell2mat(T);

         [ I N ] = size(X);           % [ 1 1500]
         [ O N ] = size(T);

        MSE00 = mean(var(t',1)) % 1.1197e+07
        MSE00a = mean(var(t',0)) %1.12041e+07

        %Normalization

        zx = zscore(cell2mat(X), 1);
        zt = zscore(cell2mat(T), 1);

            Ntrn = N-2*round(0.15*N)  
            trnind = 1:Ntrn
            Ttrn = T(trnind)

        Neq     = prod(size(Ttrn))     % 1500

        %Significant Lags were determined using the code at :

        <http://www.mathworks.com/matlabcentral/newsreader/view_thread/341287#935393>

        %sigilag95: [0:517 520 521 524 525 636:1255 1345:1384] %Significant Input Lag

        %sigflag95: [0:348 411:1207 1321:1401]  %significant Feedback lag

        rng('default')

        % %  
        % %  

FD   = 1:20; %Random Selection of sigflag subset

ID   = 1:20; %Random selection of sigilag subset

   NFD  = length(FD)   % 

   NID  = length(ID)  %

MXFD  = max(FD)      

   MXID = max(ID)

Ntrneq = prod(size(t))

   %  Nw =  ( NID*I + NFD*O + 1)*H + ( H + 1)*O

   Hub     =  -1+ceil( (Ntrneq-O) / ((NID*I)+(NFD*O)+1))  

    Hmax    =  floor(Hub/10) %  

    Hmax = 2 ==>  Nseq >>Nw :
           Hmin    = 0
           dH      = 1
           Ntrials = 25
           j=0
           rng(4151941)
           for h = Hmin:dH:Hmax
              j = j+1
              if h == 0
                  net = narxnet( ID, FD, [] );
                  Nw =  ( NID*I + NFD*O + 1)*O
              else
                  net = narxnet( ID, FD, h );
                  Nw =  ( NID*I + NFD*O + 1)*h + ( h + 1)*O
              end
              Ndof            = Ntrn-Nw
              [ Xs Xi Ai Ts ] = preparets(net,X,{},T);
              ts              = cell2mat(Ts);
              xs              = cell2mat(Xs);
              MSE00s          = mean(var(ts',1))
              MSE00as         = mean(var(ts'))
              MSEgoal         = 0.01*Ndof*MSE00as/Neq
              MinGrad         = MSEgoal/10
              net.trainParam.goal      =  MSEgoal;
              net.trainParam.min_grad  =  MinGrad;
              net.divideFcn            =  'dividetrain';
              for i = 1:Ntrials
                  net            =  configure(net,Xs,Ts);
                  [ net tr Ys ]  =  train(net,Xs,Ts,Xi,Ai);
                  ys             =  cell2mat(Ys);
                  stopcrit{i,j}  = tr.stop;
                  bestepoch(i,j) = tr.best_epoch;
                  MSE            = mse(ts-ys);
                  MSEa           = Neq*MSE/Ndof;
                  R2(i,j)        = 1-MSE/MSE00s; 
                  R2a(i,j)       = 1-MSEa/MSE00as;
              end
           end
           stopcrit   =  stopcrit    %Min grad reached (for all).
           bestepoch  =  bestepoch
           R2         =  R2
           R2a        =  R2a
           Totaltime  =  toc

[Greg Heath]

Run your code on a MATLAB dataset

 help nndatasets
 doc nndatasets

and I will compare the results with those from my code.

[Olumide Oladoyin]

I ran my code with the simplenarx_dataset

load simplenarx_dataset
X = simplenarxInputs;
T = simplenarxTargets;

and delays:

FD = 1:2;
ID = 1:2;

These were my results:

MSE00 =

      0.1021

MSE00a =

      0.1031

NFD =

       2

NID =

       2

MXFD =

       2

MXID =

       2

Hub =

      19

Hmax =

       1

Hmin =

       0

dH =

       1

Ntrials =

      10

j =

       0

j =

       1

Nw =

       5

Ndof =

      65

MSE00s =

      0.0992

MSE00as =

      0.1002

MSEgoal =

     9.3021e-04

MinGrad =

     9.3021e-05

j =

       2

Nw =

       7

Ndof =

      63

MSE00s =

      0.0992

MSE00as =

      0.1002

MSEgoal =

     9.0159e-04

MinGrad =

     9.0159e-05

stopcrit = 

      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'
      'Performance goal met.'    'Performance goal met.'

bestepoch =

       1     6
       1     4
       1     4
       1     5
       1     5
       1     7
       1     4
       1     7
       1     9
       1     7

R2 =

      0.9993    0.9983
      0.9999    0.9967
      0.9994    0.9914
      0.9999    0.9978
      0.9995    0.9980
      0.9997    0.9980
      0.9999    0.9964
      0.9996    0.9927
      1.0000    0.9909
      0.9995    0.9981

R2a =

      0.9993    0.9981
      0.9999    0.9963
      0.9993    0.9905
      0.9999    0.9975
      0.9994    0.9978
      0.9997    0.9978
      0.9998    0.9961
      0.9996    0.9919
      1.0000    0.9900
      0.9995    0.9979

Any similarity?

Thanks.

Answer 1

[Greg Heath]

I get the same results as you.

1. However, there are some code inconsistencies including:

• a. Using Ntrn = 70 with 'dividetrain'
• b. Not taking lags into account when counting equations... the equation count should be based on ttrns.

2. I now find using o (instead of s) and c as subscripts for openloop and closeloop to be more natural.

3. Sorry for not suggesting

• a. less trivial dataset examples like simpleseries and pollution.
• b. Using ID = 0

3. Remember, DIVIDETRAIN is only recommended for estimating the minimum values of ID, FD, and H that will yield acceptable performance. To get reliable estimates of performance on unseen data use data division with Ntrn as small as possible.

Hope this helps.

Greg

[Olumide Oladoyin]

Thanks Greg.