As predicted delayed outputs settle in NarX

I applied my code to data simplenarx_dataset. To do this I performed the following steps:

1 - I have done autocorrelation and cross correlation peaks to see that gives us more information. ID = 1, FD = 1

2 - I have found H, where H = 5

3 - I have created the network and have evaluated the details. Although the purpose of this post is not to evaluate the details but understand why you see a delayed response when performing closeloop, but public details and code In case of emergency there is some other error: My code is as follows(I used 80 data for training the network and 20 to check with closeloop):

p=p';

t=t';

p1=p(1:1,1:80);

p2=p(1:1,81:end);

t1=t(1,1:80);

t2=t(1,81:end);

inputSeries = tonndata(p1,true,false);

targetSeries = tonndata(t1,true,false);

inputDelays = 1:1;

feedbackDelays = 1:1;

hiddenLayerSize = 5;

net = narxnet(inputDelays,feedbackDelays,hiddenLayerSize);

[inputs,inputStates,layerStates,targets] = preparets(net,inputSeries,{},targetSeries);

net.divideFcn='divideblock';

net.divideParam.trainRatio=0.70;

net.divideParam.valRatio=0.15;

net.divideParam.testRatio=0.15;

[I N]=size(p1);

[O N]=size(t1);

N=N-1;

Neq=N*O;

ID=1;

FD=1;

Nw = (ID*I+FD*O+1)*hiddenLayerSize+(hiddenLayerSize+1)*O;

Ntrneq = N -2*round(0.15*N);

Ndof=Ntrneq-Nw;

ttotal=t1(1,1:N);

MSE00=mean(var(ttotal,1));

MSE00a=mean(var(ttotal,0));

t3=t(1,1:N);

[trainInd,valInd,testInd] = divideblock(t3,0.7,0.15,0.15);

MSEtrn00=mean(var(trainInd,1));

MSEtrn00a=mean(var(trainInd,0));

MSEval00=mean(var(valInd,1));

MSEtst00=mean(var(testInd,1));

net.trainParam.goal = 0.01*Ndof*MSEtrn00a/Ntrneq;

[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;

MSEtrn=tr.perf(end);

MSEval=tr.vperf(end);

MSEtst=tr.tperf(end);

R2trn=1-MSEtrn/MSEtrn00;

R2trna=1-MSEtrn/MSEtrn00a;

R2val=1-MSEval/MSEval00;

R2tst=1-MSEtst/MSEtst00;

and my results are:

ID=1

FD=1

H=5

N=79

Ndof=34

Neq=79

Ntrneq=55

Nw=21

O=1

I=1

R2=0.8036

R2a=0.7347

R2trn=0.8763

R2trna=0.8786

R2val=0.7862

R2tst=0.7541

As I mentioned earlier, I will not focus much on the accuracy in the answer but later will. The code I applied for closeloop was:

netc = closeloop(net);

netc.name = [net.name ' - Closed Loop'];

view(netc)

NumberOfPredictions = 15;

s=cell2mat(inputSeries);

t4=cell2mat(targetSeries);

a=s(1:1,79:80);

b=p2(1:1,1:15);

newInputSeries=[a b];

c=t4(1,80);

d=nan(1,16);

newTargetSet=[c d];

newInputSeries=tonndata(newInputSeries,true,false);

newTargetSet=tonndata(newTargetSet,true,false);

[xc,xic,aic,tc] = preparets(netc,newInputSeries,{},newTargetSet);

yPredicted = sim(netc,xc,xic,aic);

w=cell2mat(yPredicted);

plot(cell2mat(yPredicted),'DisplayName','cell2mat(yPredicted)','YdataS

ource','cell2mat(yPredicted)');figure(gcf)

plot(t2,'r','DisplayName','targetsComprobacion')

hold on

plot(w,'b','DisplayName','salidasIteradas')

title({'ITERACCIONES'})

legend('show')

hold off

and the result was the chart that you have indicated the link below where you will see it:

http://www.subirimagenes.com/otros-simenarx-8376264.html

In this picture we see the blue line (line outputs predicted) lags behind the red line (real targets). I'd like to know how I can do to get that blue line is in front of the red line, that is one step get out early. As I said, in this post I want to focus on why this happens and how I can fix it.

thank you very much

Answer 1

[Greg Heath]

% 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

[FRANCISCO]

I will describe step by step process that I followed and I will put code also.

1 - autocorrelation and partial correlation.

The code I used is as follows:

N = length(t);

ZT = zscore(t);

autocorrT = nncorr(ZT,ZT,N-1);

figure(2)

subplot(2,1,1)

autocorr(autocorrT)

title('ACF')

subplot(2,1,2)

parcorr(autocorrT)

title('PACF')

The graph obtained was presented in the following link:

http://s2.subirimagenes.com/otros/previo/thump_8382996autocorrelation.jpg

I see that the maximum lag spikes occur in 0,1,2. In this case as the Lag 2 has a higher peak, choose it as late. FD=1:2

2- crosscorrelation

The code applied is as follows:

X=p;

T=t;

[ I N ] = size(X);

[ O N ] = size(T);

Neq = N*O;

crosscorrXT0 = N*nncorr(zscore(X),zscore(T),N-1,'biased')/(N-1);

crosscorrTX0 = N*nncorr(zscore(T),zscore(X),N-1, 'biased')/(N-1);

crosscorrXT = [ crosscorrTX0(1:N-1), crosscorrXT0(N:end)];

LAG = -(N-1):(N-1);

[ crosscorrmax indmax ] = max(abs(crosscorrXT));

sigind = find(abs(crosscorrXT));

ID = LAG(sigind);

LID = length(ID(ID>0));

figure(4)

subplot(4,1,1)

plot( LAG, crosscorrXT, 'LineWidth', 2 )

title('plot')

subplot(4,1,2)

stem(LAG,crosscorrXT,'DisplayName','crosscorrXT');figure(gcf)

title('stem')

subplot(4,1,3)

plot( LAG, crosscorrXT, 'LineWidth', 2 )

xlim([0,N])

title('plot positivo')

subplot(4,1,4)

stem(LAG,crosscorrXT,'DisplayName','crosscorrXT');figure(gcf)

xlim([0,N])

title('stem positivo')

The graphs produced are:

http://s2.subirimagenes.com/otros/previo/thump_8383006crosscorrelation.jpg

I see the lag 1 is the peak that gives me more information, so ID = 1:1

3-Determination H

For the determination of hidden layers use the following code:

[I N]=size(P);

[O N]=size(T);

Neq=N*O;

Hub=floor((N-1)*O/(I+O+1)); %max H for Neq>=Nw

Hmax=floor(Hub/10);

ID = 1:1

FD = 1:2

NID= length(ID)

NFD=length(FD)

Ntrials=5;

for j=3:14

    H=j

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

    for i=1:Ntrials

        ID=1:1;

        FD=1:2;

        hiddenLayerSize=H;

        rng(0)

        net=narxnet(ID,FD,hiddenLayerSize);
        [inputs,inputstates,layerstates,targets]=preparets(net,P,{},T);

        [ I N ] = size(inputs)

        [ O N ] = size(targets)

        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

        ttrn=cell2mat(ttrn);

        ytrn00= mean(ttrn,2);  

        Nw00 = size(ytrn00,2)

        Ndof00 = Ntrneq-Nw00

        MSEtrn00=mean(var(ttrn,1))

        MSEtrn00a=mean(var(ttrn,0))

        tval=cell2mat(tval);

        MSEval00 = mse(tval-ytrn00)

        ttest=cell2mat(ttest);

        MSEtst00 = mse(ttest-ytrn00)

        net.trainParam.goal = 0.01*Ndof*MSEtrn00a/Ntrneq;

        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)

        MSEtrn=tr.perf(end)

        MSEtrna = Ntrneq*MSEtrn/Ndof

        MSEval=tr.vperf(end)

        MSEtst=tr.tperf(end)

        R2trn=1-MSEtrn/MSEtrn00

        R2trna=1-MSEtrna/MSEtrn00a

        R2val=1-MSEval/MSEval00

        R2tst=1-MSEtst/MSEtst00

        precisiones{i,j}=[R2trn R2trna R2val R2tst];

    end
end

I note that the best accuracies are obtained with H = 7;

R2trn = 0.9993

R2trna = 0.9981

R2val = 0.9992;

R2tst = 0.9976

4-Creation of the network

I use the following code:

[P, T ] = simplenarx_dataset;

whos

p= cell2mat(P);

t = cell2mat(T);

ID = 1:1

FD = 1:2

H = 7

rng(0)

net = narxnet(ID,FD,H);

[inputs,inputStates,layerStates,targets] = preparets(net,P,{},T);

whos P T inputs inputStates layerStates targets

[ I N ] = size(inputs)

[ O N ] = size(targets)

NID= length(ID)

NFD=length(FD)

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

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

ttrn=cell2mat(ttrn);

ytrn00= mean(ttrn,2);

Nw00 = size(ytrn00,2)

Ndof00 = Ntrneq-Nw00

MSEtrn00=mean(var(ttrn,1))

MSEtrn00a=mean(var(ttrn,0))

tval=cell2mat(tval);

MSEval00 = mse(tval-ytrn00)

ttest=cell2mat(ttest);

MSEtst00 = mse(ttest-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);

MSEtrn=tr.perf(end)

MSEtrna = Ntrneq*MSEtrn/Ndof

MSEval=tr.vperf(end)

MSEtst=tr.tperf(end)

R2trn=1-MSEtrn/MSEtrn00

R2trna=1-MSEtrna/MSEtrn00a

R2val=1-MSEval/MSEval00

R2tst=1-MSEtst/MSEtst00

plot(cell2mat(outputs),'b','DisplayName','outputs')

hold on

plot(cell2mat(T),'r','DisplayName','targetSeries')

title({'ENTRENAMIENTO'})

legend('show')

hold off

The training is graphic:

http://s2.subirimagenes.com/otros/previo/thump_8383099creationnetwork.jpg

and the results obtained are:

N=98

NFD=2

NID=1

Ndof=27

Ndof00=69

Ntrneq=70

Nw=43

R2trn=0.9993

R2trna=0.9981

R2val=0.9992

R2tst=0.9976

The outputs in closeloop still appear delayed me, but I wish you gave me your opinion about this that I've done now.

thank you very much

Answer 2

[Greg Heath]

    1. I do not have autocorr or parcorr AND your plots are too fuzzy
    for me to understand.

    2.Since training parameters should only depend on training data,
    I use ztrn =zscore(ttrn,1), autocorrt = nncorr(ztrn,ztrn,Ntrn-1,'biased')
    and find a 95% confidence level for significant auto AND
    cross-correlations with abs values  >= 0.14 by repeating 100
    times and averaging:
    a. Crosscorrelate ztrn with n = zscore(randn(1,Ntrn),1)
    b. Correct the nncorr symmetry bug by concatenating: 
       crosscorrzn = [ crosscorrnz(1:Ntrn-1) crosscorrzn(Ntrn:end)]
    c. Sort the absolute values and find the significance threshold
    as the value at index floor( 0.95*(2*Ntrn-1))

    3. You used crosscorrxt instead of crosscorrtx. You found the
    max of crosscorrxt instead of all abs values >= the significance
    threshold for nonnegative lags (indices >= Ntrn) of crosscorrtx.

    THIS MAY ACCOUNT FOR YOUR ERRONEOUS SHIFT

    4. The ideal lags are ID = 1, FD = 1,2

    5. To determine Hub use Ntrneq, NID and NFD (see your equation
    for Nw)

        6. I typically find Ntrials = 5 to be too small

        7. You are using divideblock inside the inner loop. Since the data
        division doesn't change, use divideblock before the outer loop
        to obtain the indices and the data splits ttrn,tval,tst, etc.

        8. Inside the inner loop just assign the indices to
        net.divide.trainInd etc.

        9. Continuously find and save the best net within the inner loop
        after you have R2val to compare with the current R2valmax.
        Then when you exit the outer loop you will already have the best
        net (and all of it's statistics) to use for the closeloop stage.

        10. Test the netc on all of the data and overlay plot yc on tsc.

        11. Train netc and compare trn/val/tst results with the openloop
        solution and the untrained closeloop solution.

        12. Unfortunately, I have found that a good closeloop solution
        depends very critically on the choices of ID, FD and H. Just
        accepting the defaults is very suboptimal.

        Hope this helps.

• Thank you for formally accepting my answer*

Greg

[FRANCISCO]

sorry for this long without commenting. I've been testing, but I think I still have the same problem unclear. The ultimate end of my network is to predict the value of the next day, today introduced inputs. How could I do this with the closed loop pump? since I already tried openloop network, I made closed loop pump with the training data, and now how would that prediction as if in a real way?. I have seen that using for example ID = 2, DF = 2, the network response is delayed two rejections, whereas if uses cookies for example ID = 1 DF = 2 the network response is ahead in time but requires me to enter futulos inputs (which is impossible). How could I make the prediction in shape?.

thank you very much

Answer 3

[Greg Heath]

1. Your use of the words pump, rejections, cookies and futulos(not English) are not appropriate and somewhat confusing.

2. Once you have an openloop solution, y, netc = closeloop(net); followed by preparets and yc = netc(Xsc,Xic,Aic) should yield your closeloop output.

3. If yc and y differ substantially, you can either continue training the openloop configuration and repeat OR, directly train netc.

4. As long as ID > 0, you have a predictive net. So, I think you are misinterpreting your output.

[Constantin Silkin]

I'm using NARX and have just the same problem. Actually the network doesn't carry out a prediction. Instead it repeats slightly changed previous value. Nevertheless, I could organize a prediction! For this purpose I prolong TargetSeries and InputSeries. I just duplicate the last value of these series. In fact I give for network own option of a prediction by the principle "tomorrow will be same as today". The network processes this naive forecast according to own knowledge and corrects it. Accuracy of such prediction is about 70%. Whether I am right?

[Greg Heath]

Sorry, I do not understand.

Greg