Auto correlation and cross correlation with NARXnet to find ID and FD

Hello,

I have recently started working on neural networks and time series with Matlab. So far, I have been using NARX networks and have been training them and then looking at the inerrcorr and errcorr plots to know if I needed to use bigger delays but it has come to my attention that there is a more "rigorous" way of doing it using the nncorr function.

I have looked at many questions and answers but there are still a few steps along the process that I don't understand. I am going to use a code derived from the one found in the topic NARXnet non-integer delays problem to show which parts I don't get:

X = inputs; % 1x2323 mat
T = targets; % 1x2323 mat

N = lenght(T)

crosscorrXT = nncorr(X,T,N-1);
autocorrT = nncorr(T,T,N-1);
crosscorrXT(1:N-1) = []; % Delete negative delays
autocorrT(1:N-1) = [];
sigthresh95 = 0.21 % Significance threshold
sigcrossind = crosscorrXT( crosscorrXT >= sigthresh95 )
sigautoind = autocorrT( autocorrT >= sigthresh95 )
inputdelays = sigcrossind(sigcrossind <= 35)
feedbackdelays = sigautoind(sigautoind <= 35)
feedbackdelays(1)=[] % Delete zero delay

Though I understand from the answer given there that the value used isn't correct, I still don't understand how to get the 95% confidence limit. I also don't understand what the last three lines of her/his code do.

I found on another topic another way of getting sigthresh95:

sigthresh95 = crosscorrXT(floor(0.95*(2*N-1)))

but this line returns the error " Attempted to access crosscorrXT(4412); index out of bounds because numel(crosscorrXT)=2323" so clearly I am not using it correctly and I don't understand how it is supposed to work.

I think I'm not not clear on what the 95% is of.

Answer 1

[Greg Heath]

Look for my recent NEWSGROUP posts (2015)

 greg nncorr 2015

Especially

 8 Jul 2015 OPENLOOP NARNET TUTORIAL EXAMPLE
 2 Jun 2015 SIGNIFICANT AUTOCORRELATION DELAYS OF THE SIMPLENAR_DATASET
 13 Apr 2015 NARNET TUTORIAL ON MULTISTEP AHEAD PREDICTIONS: MESSAGES 1, 7, 11

Hope this helps.

Thank you for formally accepting my answer

Greg

[Nikita OKS]

Dear Greg,

Thank you for your help and for the suggested posts. I found the first two right after posting my question here and I used them to write my own code but here is what it returns:

[I N] = size(inputs); % [1 2323]
[O N] = size(targets); % [1 2323]
MSE00 = mean(var(targets',1)) % 8.0153e+08
Ntrn = N-2*round(0.15*N) % 1627
Ntrneq = Ntrn*O % 1627
H = 25;
Nw = (I+1)*H+(H+1)*O % 76
r = Ntrneq/Nw % 21.4079
MSEgoal = 0.01*MSE00 % 8.0153e+06

and

[O N] = size(targets); % [1 2323]
zt = zscore(targets',1)';
L = floor (0.95*(2*N-1)); % 4412

rng('default')

for i = 1:100

      zn=zscore(randn(1,N));
      autocorrn = nncorr(zn,zn,N-1,'biased');
      sortabsautocorrn = sort(abs(autocorrn));
      thresh95(i) = sortabsautocorrn(L);
  end

meanthresh95 = mean(thresh95); %0.0304

autocorrt = nncorr(zt,zt, N-1,'biased');
sigthresh = find(abs(autocorrt(N+1:2*N-1))>=meanthresh95)
% sigthresh 1x1953 mat
% sigthresh = [1:660 742:2035]

FD = sigthresh;
NFD = numel(FD); % 1953
Ntrn = N-2*round(0.15*N) % 1627
Ntrneq = Ntrn*O % 1627
Hub = floor((Ntrneq-O)/(NFD*O+O+1)) % 0

As you can tell, my sigthresh is very large and I don't see how I could use it for my delays ( and I haven't even tried to create a double-loop net yet because of the Hub equal to zero).

Because I am not sure to understand some parts of this code (especially how the hub is calculated or how the for loop is used to compute the confidence limit), I fail to write a code that works with my own data.

Finally, I also noticed a slight difference in the way you obtained sigthresh in the two posts as one defines it as shown above when the other adds a -1 to the formula and I am not sure to understand why.

sigthresh = -1+find(abs(autocorrt(N+1:2*N-1))>=meanthresh95)

Again, thank you for your help

Nikita

Answer 2

[Greg Heath]

 [ I N ] = size(inputs); % [1 2323]
 [ O N ] = size(targets); % [1 2323]
 MSE00   = mean(var(targets',1)) % 8.0153e+08

%GEH1: With a value this huge, why didn't you normalize??

 Ntrn   = N-2*round(0.15*N) % 1627
 Ntrneq = Ntrn*O % 1627
 H      = 25;

%GEH2: Why H so high?? Typically, you want H and Nw to be as small as possible to increase the stability, robustness and generalization capability of the design. At this point I would suggest starting with the default H=10.

 Nw      = (I+1)*H+(H+1)*O % 76
 r       = Ntrneq/Nw          % 21.408
 MSEgoal = 0.01*MSE00       % 8015300

% GEH3: This only leads to an Rsquared of 0.990. For an openloop net, It is probably better to try for 0.995 or 0.999

 [O N ] = size(targets);        % [1 2323]
 zt     = zscore(targets',1)';
 L      = floor (0.95*(2*N-1)); % 4412

% GEH4: To keep validation and test subset predictions as unbiased as possible, probably should only use the training or, if needed, design (training+validation) subsets to determine significant lags and number of hidden nodes. e.g., use Ntrn above.

 rng('default')
 for i = 1:100

      zn = zscore(randn(1,N));

% GEH5: Use biased form : zscore(..., 1)

   autocorrn        = nncorr(zn,zn,N-1,'biased');
   sortabsautocorrn = sort(abs(autocorrn));
   thresh95(i)      = sortabsautocorrn(L);
  end
meanthresh95 = mean(thresh95); %0.0304

 autocorrt = nncorr(zt,zt, N-1,'biased');
 sigthresh = find(abs(autocorrt(N+1:2*N-1))>=meanthresh95)

% sigthresh 1x1953 mat % sigthresh = [1:660 742:2035]

% GEH6: Should use Ntrn and ttrn in all above calculations

% GEH7: For Narxnet also need crosscorrxt calculations.

 FD = sigthresh;
 NFD = numel(FD); % 1953

% GEH8: No, No, No! Only use a subset of FD long enough to achieve MSEgoal.

 Ntrn   = N-2*round(0.15*N) % 1627
 Ntrneq = Ntrn*O % 1627
 Hub    = floor((Ntrneq-O)/(NFD*O+O+1)) % 0

As you can tell, my sigthresh is very large and I don't see how I could use it for my delays ( and I haven't even tried to create a double-loop net yet because of the Hub equal to zero).

% See GEH8

Because I am not sure to understand some parts of this code (especially how the hub is calculated or how the for loop is used to compute the confidence limit), I fail to write a code that works with my own data.

%GEH9: Hub is the upper bound for Nw <= Ntrneq (No. of unknowns no greater than the number of training equations)

% GEH10: See your stats refererence or wikipedia re defining statistical significance w.r.t. noise statistics ( Remember p- values etc?)

Finally, I also noticed a slight difference in the way you obtained sigthresh in the two posts as one defines it as shown above when the other adds a -1 to the formula and I am not sure to understand why. sigthresh = -1+find(abs(autocorrt(N+1:2*N-1))>=meanthresh95)

% GEH11: Cannot have FD = 0 when loop is closed. But can have ID=0

 sigthreshxcorr = -1+ find(abs(crosscorrxt(N:2*N-1))>=meanthresh95)
 sigthreshacorr = find(abs(autocorrt(N+1:2*N-1))>=meanthresh95)

Hope this helps.

Greg

[Nikita OKS]

Dear Greg,

I have modified my code and it does yield better results now:

% tic
[norm_inputs,xs] = mapminmax(inputs);
[norm_targets,ts] = mapminmax(targets);

[O N] = size(norm_targets); % [1 2323]
[I N] = size(norm_inputs);
zt = zscore(norm_targets',1)';
zx = zscore(norm_inputs',1)';

MSE00 = mean(var(t',1)) % 0.1727
MSEgoal = 0.005*MSE00 % 0.0017
Ntrn = N-2*round(0.15*N) % 1627
Ntrneq = Ntrn*O % 1627
L = floor (0.95*(2*Ntrn-1)); % 3090

rng('default')

for i = 1:100

      zn=zscore(randn(1,Ntrn),1);
      autocorrn = nncorr(zn,zn,Ntrn-1,'biased');
      sortabsautocorrn = sort(abs(autocorrn));
      thresh95(i) = sortabsautocorrn(L);
  end

meanthresh95 = mean(thresh95); %0.0362

autocorrt = nncorr(zt,zt, Ntrn-1,'biased');
sigflag95 = find(abs(autocorrt(Ntrn+1:2*Ntrn-1))>=meanthresh95);
% sigflag95 1x1540 mat
% sigflag95 = [1:659 746:1626]

crosscorrxt = nncorr(zx,zt, Ntrn-1,'biased');
sigilag95 = -1 + find(abs(crosscorrxt(Ntrn+1:2*Ntrn-1))>=meanthresh95);
% sigilag95 1x1540 mat
% sigilag95 = [1:657 744:1625]

FD = 1:10; % ??? See my question below
NFD = numel(FD); % 10

ID = 1:10;
NID = numel(ID);

% Nw = (NID*I+NFD*O+1)*H+(H+1)*O % No. of unknown weights
% Ntrneq >= Nw <==> H <= Hub % Unknowns don't outnumber equations
Hub = floor((Ntrneq-O)/(NID*I+NFD*O+O+1)) % 73
% I know it is only an upper bound but it still seems extreme

At this point, I tried to find a good subset of FD/ID using the following network built in a loop:

X = tonndata(train_inputs,true,false);
% train_inputs 1x1627 mat with 1627 first values from norm_inputs
T = tonndata(train_targets,true,false);
% train_targets 1x1627 mat with 1627 first values from norm_targets

trainFcn = 'trainbr';

for j = 1:1:5
i= 2
performance = 1;

      while performance > 0.0017

          inputDelays = 1:i;
          feedbackDelays = 1:i;
          hiddenLayerSize = 8; % Hub
          net = narxnet(inputDelays,feedbackDelays,hiddenLayerSize,'open',trainFcn);

          net.inputs{1}.processFcns = {'removeconstantrows', 'mapminmax'};
          net.inputs{2}.processFcns = {'removeconstantrows', 'mapminmax'};

          [x,xi,ai,t] = preparets(net,X,{},T);

          net.divideFcn = 'dividetrain';

          net.performFcn = 'mse';
          net.trainParam.epochs = 5000;

          net.plotFcns = {'plotperform','plottrainstate', 'ploterrhist', ...
              'plotregression', 'plotresponse', 'ploterrcorr', 'plotinerrcorr'};

          [net,tr] = train(net,x,t,xi,ai);

          y = net(x,xi,ai);
          e = gsubtract(t,y);
          performance = perform(net,t,y)

          i=i+1
      end
  end

Now my problem is that I never achieve the 0.0017 for the same delay. Of all the five networks created within the for loop, each achieved MSEgoal for FDs and IDs of 10, 8, 6, 8, 9 respectively and I am not quite sure of which one I should go with (the smallest and have my final network train enough times to have it achieve the desired performance ?).

Finally, here is the rest of the code where I try to find the optimal size for my hidden layer (I haven't run it yet as I am still working on the delays).

Hmax = 8;
dH = 1;
Hmin = 0;
Ntrials = 10;

X = tonndata(train_inputs,true,false); 
train_inputs 1x1627 mat with 1627 first values from norm_inputs
T = tonndata(train_targets,true,false);
train_targets 1x1627 mat with 1627 first values from norm_targets

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,'open',trainFcn);
        Nw = (NID*I+NFD*O+1)*h+(h+1)*O;
    end

      Ndof = Ntrneq - Nw % <= 0 for h >= 8
      neto.divideFcn = 'dividetrain';
      neto.performFcn = 'mse';

      neto.inputs{1}.processFcns = {'removeconstantrows', 'mapminmax'};
      neto.inputs{2}.processFcns = {'removeconstantrows', 'mapminmax'};
      [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] = train(neto,Xo,To,Xoi,Aoi);

          % Yo = net(Xo,Xoi,Aoi);
          % Eo = gsubtract(To,Yo);

          % performance = perform(neto,To,Yo);
          stopcrit{i,j} = tro.stop;
          R2o(i,j) = 1 - mse(Eo)/MSE00o;
      end
  end

result = [(Hmin:h:Hmax); R2o]
stopcrit = stopcrit
elapsedtime = toc

Sorry about this long post which hopefully will be my last.

Thank you,

Nikita

[Greg Heath]

Don't be sorry !! I enjoyed it because it reminded me of a few things that I forgot to include in some tutorials AND the joint optimization problem with partial correlations which I had sidestepped.

There are still a few things about trainbr that bother me so, ... be careful.

 1. Begin with zscore transformation of training data 
 2. Perform Outlier & BadData checks using zscore parameters: 
    a. Apply (max-trnmean)/trnstd & (trnmean-min)/trnstd to trn/val/tst
    b. Plots 
 3. You can use ID >= 0 but FD > 0
 4. Use Ntrn to find ID and FD. Only if Ntrn is small should you  ever consider using Ntrn+Nval. Otherwise you lose the credibility of assuming val performance estimates are relatively unbiased.
 5. A huge Hub just means you have a lot of data and probably don't have to worry too much about overtraining an overfit net. 
 6. a. I have never looped to optimize ID and FD. Instead, I plot autocorr and crosscorr in blue or black with the significant correlations points overlayed in red circles 'ro'. 
   b. Sometimes looking at the plots have caused me to use the median or mean+std to define significance.
 7. a.Why did you feel the following  statements were necessary?
     trainFcn = 'trainbr'
     net.inputs{1}.processFcns = {'removeconstantrows', 'mapminmax'};
     net.inputs{2}.processFcns = {'removeconstantrows', 'mapminmax'};
     net.divideFcn = 'dividetrain';
     net.performFcn = 'mse';
     net.trainParam.epochs = 5000;
net.plotFcns = {'plotperform','plottrainstate', 'ploterrhist', ...
              'plotregression', 'plotresponse', 'ploterrcorr', 'plotinerrcorr'};
   b. Are trainbr and dividetrain compatible?
 8. You have spent quite a bit of time trying to optimize delays given a fixed H. Then planning to fix those delays and optimizing H. Interesting because I had decided to go the other way. BUT,  after I have gotten a reasonable combination, I just accepted it and didn't worry about a joint optimum solution. 
 9. Be aware that there is a flaw in this type of joint optimality search: The straightforward correlation functions do not take interactions into account. For example, assume the top two significant delays are 1 and 4, respectively.  If the 1 delay correlation is removed by subtraction, the remaining top significant correlation delay may no longer be 4.

 This is taken into account using the the theory of PARTIAL CORRELATIONS, the theory of which I am not familiar. However, it probably is not much different than the familiar stepwise and stagewise methods for selection of regression variables. 

 I think the Econometrics Tool Box does use partial correlations for stepwise or stagewise variable selection.
 10. You did not include the final delay states xf and af in your outputs.

Hope this helps.

Thank you for formally accepting my answer

Greg

[Nikita OKS]

Dear Greg,

Thank you for your reply.

   b. Are trainbr and dividetrain compatible?

They seem to be, yes. But I do have a question about dividetrain as I remember seeing in the user's manual that they left the divideFcn parameter blank (i.e. just the quotes) and I don't know if there's a difference between not telling the code how to divide the series and telling it to just use all the data for training.

Anyway, I haven't had time to add the modifications you suggested to my code yet as something else has required my attention but I will as soon as I can and I think I should be good after that.

Thank you for your much appreciated help.

Nikita

Answer 3

[mohamad Macabangin]

I am new here. I have been searching for greg's NEWSGROUP posts (2015) and I always ended up with Matlab answers post. Can I have the link?