MATLAB: Finding Optimal ID, FD and Hidden Nodes for NARXNET

hidden nodeinput delaynarxnetneural networkoutput delay

Hi, After a lengthy research, I have finally have better understanding about ID and FD. I have then put some code together to find the optimal ID and FD and then using these ID and FD to find optimal hidden node for my NARXNet using simplenarx dataset. While I am convince that it is correct but I am not very confident if this is the correct way of doing things so I would really appreciate any comments/correction if any.

Some additional question are: 1) Should I use data division such as 60/20/20 in the double for loop? 2) I used intersect command to find the subset of lags, is this correct way of doing it?

Below are my codes:

if true
  % code
end
close all, clear all, clc, 
tic
plt=0;
[X,T] = simplenarx_dataset; % simplenarx_dataset;
x = cell2mat(X);
t = cell2mat(T);
[ I N ] = size(x); % [ 1 100 ]

[ O N ] = size(t); % [ 1 100 ]
MSE00 = var(t',1) % 0.1021
% Calculate Z-Score for input (x) and target (t)
zx = zscore(x, 1);
zt = zscore(t, 1);
% Plot Input & Output for both original and transformed (Z-scored)
plt = plt+1,figure(plt);
subplot(221)
plot(x)
title('SIMPLENARX INPUT SERIES')
subplot(222)
plot(zx)
title('STANDARDIZED INPUT SERIES')
subplot(223)
plot(t)
title('SIMPLENARX OUTPUT SERIES')
subplot(224)
plot(zt)
title('STANDARDIZED OUTPUT SERIES')
rng('default')
L = floor(0.95*(2*N-1)) % 189
for i = 1:100 % 1: length of the data
  n = zscore(randn(1,N),1);
  autocorrn = nncorr( n,n, N-1, 'biased');
  sortabsautocorrn = sort(abs(autocorrn));
  thresh95(i) = sortabsautocorrn(L);
end
sigthresh95 = mean(thresh95) % 0.1517
minthresh95 = min(thresh95) % 0.1139 
medthresh95 = median(thresh95) % 0.1497
stdthresh95 = std(thresh95) % 0.0234
maxthresh95 = max(thresh95) % 0.2321 
%%CORRELATIONS
%%%%%TARGET AUTOCORRELATION %%%%%%%
% 



autocorrt = nncorr(zt,zt,N-1,'biased');
sigflag95 = -1+ find(abs(autocorrt(N:2*N-1))>=sigthresh95) %significant Feedback Delay (FD) => [0 2 3 4 5 7 9 10 12 67 69]
% 
plt = plt+1, figure(plt);
hold on
plot(0:N-1, -sigthresh95*ones(1,N),'b--')
plot(0:N-1, zeros(1,N),'k')
plot(0:N-1, sigthresh95*ones(1,N),'b--')
plot(0:N-1, autocorrt(N:2*N-1))
plot(sigflag95,autocorrt(N+sigflag95),'ro')
title('SIGNIFICANT TARGET AUTOCORRELATIONS (FD)')
%





%%%%%%INPUT-TARGET CROSSCORRELATION %%%%%%
%
crosscorrxt = nncorr(zx,zt,N-1,'biased');
sigilag95 = -1 + find(abs(crosscorrxt(N:2*N-1))>=sigthresh95) %significant Input Delay (ID) => [0 1 3 4 5 6 8 10 13 17]
% 
plt = plt+1, figure(plt);
hold on
plot(0:N-1, -sigthresh95*ones(1,N),'b--')
plot(0:N-1, zeros(1,N),'k')
plot(0:N-1, sigthresh95*ones(1,N),'b--')
plot(0:N-1, crosscorrxt(N:2*N-1))
plot(sigilag95,crosscorrxt(N+sigilag95),'ro')
title('SIGNIFICANT INPUT-TARGET CROSSCORRELATIONS (ID)')
%%Using Fixed ID and FD to Find Optimal Number of Hidden Node
%
subset_ID_FD = intersect(sigflag95, sigilag95)
Opti_ID_FD = max(subset_ID_FD);
Ntrn = N-2*round(0.15*N) % default 0.7/0.15/0.15 trn/val/tst ratios
trnind = 1:Ntrn;
Ttrn = T(trnind);
Ntrneq = prod(size(Ttrn)) % Product of element
%ID = 1:2 %default for Prediction
ID = 1:Opti_ID_FD; % 0:2 % Regression (default)
FD = 1:Opti_ID_FD; % 1:2 % default (default)
NID = length(ID); % 10

NFD = length(FD); % 10
LDB = max([ID,FD]) % Length of the delay buffer = 10
Hub = floor((Ntrneq-O)/(NFD*O+O+1)) % 5
Hmax = Hub; % 2 is sufficient to get R2=0.999
dH =1;
Hmin = 1;
Ntrials = 10;
%
trainFcn = 'trainbr'
%
rng('default')
j=0
for h = Hmin:dH:Hmax
  j=j+1
  if h==0
      neto = narxnet(ID,FD,[],'open',trainFcn);
      Nw = (NID*I+NFD*O+1)*h+(h+1)*O;
  else
      neto = narxnet(ID,FD,h);
      Nw = (NID*I+NFD*O+1)*h+(h+1)*O;
  end
  Ndof = Ntrneq-Nw % Ndof <=0 for H >= 4
  neto.divideFcn = 'dividetrain'; % No data division
  neto.performFcn = 'mse';
  [Xo Xoi Aoi To ] = preparets(neto,X,{},T);
  to = cell2mat(To);
  MSE00o = var(to,1)
  MSE00oa = var(to,0)
  MSEgoal = 0.005*max(Ndof,0)*MSE00oa/Ntrneq
  MinGrad = MSEgoal/100
  neto.trainParam.goal = MSEgoal;
  neto.trainParam.min_grad = MinGrad;
%   
  for i= 1:Ntrials
      % Save state of RNG for duplication
      s(i) = rng; 
      neto = configure(neto,Xo,To);
      [neto tro Yo Eo Xof Aof ] = train(neto,Xo, To, Xoi, Aoi);
      % Eo = gsubtract(To,Yo);
      % stopcrit{i,j} = tro.stop;
      R2o(i,j) = 1 - mse(Eo)/MSE00o;
  end
end
result = [ (Hmin:dH:Hmax); R2o ]
% stopcrit = stopcrit;
elapsedtime = toc % 194.8385
%
% result =
% 
%     1.0000    2.0000    3.0000    4.0000    5.0000
%     0.9966    0.9996    0.9999    1.0000    1.0000
%     0.9968    0.9990    0.9998    1.0000    1.0000
%     0.9967    0.9990    0.9998    1.0000    1.0000
%     0.9967    0.9991    0.9998    1.0000    1.0000
%     0.9968    0.9999    0.9999    1.0000    1.0000
%     0.9974    0.9990    1.0000    1.0000    1.0000
%     0.9970    0.9985    0.9998    1.0000    1.0000
%     0.9976    0.9987    0.9999    1.0000    1.0000
%     0.9986    0.9990    0.9998    1.0000    1.0000
%     0.9967    0.9986    0.9999    1.0000    1.0000

Many thanks! Teck

Best Answer

> Hi, After a lengthy research, I have finally have better Understanding about ID and FD. I have then put some code together to find the optimal ID and FD and then using these ID and FD to find optimal hidden node for my NARXNet using simplenarx dataset. While I am convince that it is correct but I am not very Confident if this is the correct way of doing things so I would really appreciate any comments/correction if any.

> Some additional question are:

1) Should I use data division such as 60/20/20 in the double for loop?

 You have used TRAINBR for which Nval = 0 and performFcn = msereg
 HOWEVER,you have imposed performFcn = mse and DIVIDETRAIN for which Ntst = 0 .
 VERY CONFUSING!

2) I used intersect command to find the subset of lags, is this correct way of doing it?

 NO.
 GEH1 'WHAT IS THE FOLLOWING COMMAND FOR?'

if true % code end

 GEH2 = 'USE ONLY TRAINING DATA TO DETERMINE DELAYS' 
 GEH3 = 'REMAINING POST STRICTLY VALID ONLY FOR I = O = 1 '
       ' MODIFICATIONS NEEDED FOR MULTIVARIABLE DATA'
 GEH4 = 'CANNOT USE FD = 0 WILL GET ERROR'

> Using Fixed ID and FD to Find Optimal Number of Hidden Node

subset_ID_FD = intersect(sigflag95, sigilag95)

GEH5 =  '0     3     4     5    10'

Opti_ID_FD = max(subset_ID_FD);

GEH6 = 'Opti_ID_FD = 10 NOT NECESSARILY OPTIMAL!'

Ntrn = N-2*round(0.15*N) % default 0.7/0.15/0.15 trn/val/tst ratios

GEH7 = 'ABOVE NOT VALID FOR TRAINBR  0.85/0/0.15'

%ID = 1:2 %default for Prediction ID = 1:Opti_ID_FD; % 0:2 % Regression (default)

 GEH8 = 'ZERO DELAY IS NOT A MATLAB DEFAULT'

Hub = floor((Ntrneq-O)/(NFD*O+O+1)) % 5

 GEH9 = 'Hub = (Ntrneq-O)/(NID*I+NFD*O+1)= 3.29'
 Hmax = Hub; % 2 is sufficient to get R2=0.999
dH =1;
Hmin = 1;
Ntrials = 10;
%

trainFcn = 'trainbr'
%
rng('default')
j=0
for h = Hmin:dH:Hmax
  j=j+1
  if h==0
      neto = narxnet(ID,FD,[],'open',trainFcn);
      Nw = (NID*I+NFD*O+1)*h+(h+1)*O;
 GEH10 = 'Nw = (NID+1)*O'

neto.divideFcn = 'dividetrain'; % No data division

 GEH11 = 'NEED NONTRAINING DATA FOR UNBIASED PREDICTION !!'

neto.performFcn = 'mse';

 GEH12 ' TRAINBR USES MSEREG WITH NO VAL SET'
 [Xo Xoi Aoi To ] = preparets(neto,X,{},T);
 to = cell2mat(To);
 MSE00o = var(to,1)
 MSE00oa = var(to,0)
 MSEgoal = 0.005*max(Ndof,0)*MSE00oa/Ntrneq
 GEH13 = 'ONLY USE TRAINING DATA TO COMPUTE TRAINING PARAMETERS'
 GEH14 = SUGGEST LOOKING AT TRAINING RECORD tro.

Hope this helps.

 *Thank you for formally accepting my answer*

Greg

Related Solutions

MATLAB: As predicted delayed outputs settle in NarX

% 1. Selected ending semicolons can be removed to aid debugging

 [P, T ] = simplenarx_dataset;
 whos
 p= cell2mat(P);
 t = cell2mat(T);
 ID = 1:1
 FD = 1:1
 H = 5
 NID= length(ID)
 NFD=length(FD)
 Nw = (NID*I+NFD*O+1)*H+(H+1)*O

% 2. Use NID and NFD for Nw in case delays are not single

% 3. No need to use tonndata because the simplenarx_data set is instantly ready for preparets.

% 4. No need for (p1,t1) and (p2,t2). Delete both.

% 5. Input delays are suboptimal. Did you try to find the significant lags of the target/input cross-correlation function?

% 6. Feedback delays are suboptimal. Did you try to find the significant lags of the target autocorrelation function?

% 7. H is suboptimal. Was it chosen using the suboptimal delays? If so, please explain how.

 rng(0)
 net = narxnet(ID,FD,H);
 [inputs,inputStates,layerStates,targets] = preparets(net,P,{},T);
 whos P T inputs inputStates layerStates targets

%8. N=N-1: DELETE. NOT A GOOD IDEA TO USE A VARIABLE OR PARAMETER NAME ON BOTH SIDES OF AN EQUATION. BESIDES, PREPARETS OUTPUTS THE CORRECT DIMENSIONS

 [ I N ] = size(inputs)
 [ O N ] = size(targets)

% 9. No need for ttotal it should be the same as targets. % No need for Neq, MSE00,MSE00a and t4. Delete

 net.divideFcn='divideblock';
 [trainInd,valInd,testInd] = divideblock(N,0.7,0.15,0.15); 
 ttrn = targets(trainInd);
 tval = targets(valInd);
 ttest = targets(testInd);
 Ntst = length(ttrn) 
 Nval = length(valInd)
 Ntst = length(testInd)
 Ntrneq = prod(size(ttrn))     % Ntrn*O
 Ndof = Ntrneq-Nw

%Naive Constant Output Model

 ytrn00= mean(ttrn,2);     
 Nw00 = size(ytrn00,2)
 Ndof00 = Ntrneq-Nw00
 MSEtrn00 = sse(ttrn-ytn000)/Ntrneq
 MSEtrn00=mean(var(ttrn,1))
 MSEtrn00a = sse(ttrn-ytrn00)/Ndof00
 MSEtrn00a=mean(var(ttrn,0))

%9. MSEval00 and MSEtst00 should be obtained from the Naive constant output model output

 MSEval00 = mse(tval-ytrn00)
 MSEtst00 = mse(tttst-ytrn00)
 net.trainParam.goal = 0.01*Ndof*MSEtrn00a/Ntrneq;  % R2trna >= 0.99
 rng(0)
 [net,tr,Ys,Es,Xf,Af] = train(net,inputs,targets,inputStates,layerStates);
 outputs = net(inputs,inputStates,layerStates);
 errors = gsubtract(targets,outputs);
 MSE = perform(net,targets,outputs);
 MSEa=Neq*MSE/(Neq-Nw)
 R2=1-MSE/MSE00
 R2a=1-MSEa/MSE00a

% 10. The DOF "a"djustment is only applied to the training data % 11. Can delete the last 6 equations that refer to all of the data instead of the trn/val/tst division.

 MSEtrn=tr.perf(end)
 MSEtrna = Ntrneq*MSEtrn/Ndof
 MSEval=tr.vperf(end)
 MSEtst=tr.tperf(end)

% 12.Using "end" is only valid if training converges because of tr.min_grad (not valstop ). Better to use "tr.best_epoch".

 R2trn=1-MSEtrn/MSEtrn00
 R2trna=1-MSEtrna/MSEtrn00a

%13 Original MSEtrna misprint.

 R2val=1-MSEval/MSEval00
 R2tst=1-MSEtst/MSEtst00

and my results are:

% 14. Unable to compare results because you did not intialize the RNG before the first call of the RNG in the net creation command net = ... where H =5 random weights were assigned to the input bias. I will use rng(0) before the net creation.

I will do the closeloop part next.

Thank you for formally accepting my answer

Greg

MATLAB: Input and feedback delay in narxnet

 1. Ntrn ~  0.7*651 ~  456 
 2. Use Ntrn for significant correlation threshold estimates
 3. 1000 repetitions seems extreme. How much difference if 
    you only  use 100 or 200?
 4. Mislabeled crosscorr as autocorr
 5. Print out all significant correlations. Do not try to 
    read them from the graphs.

Hope this helps.

Thank you for formally accepting my answer*

Greg

Best Answer

Related Solutions

MATLAB: As predicted delayed outputs settle in NarX

MATLAB: Input and feedback delay in narxnet

Related Question