

clear all
close all
randn('state',sum(100*clock))

niter=300;
pRKHS=40; % length of eigenfunctions
p=20;%input('# of predictor i.r. samples = ');%lunghezza risposte impulsive predittore (100)
p2=300;%numero dati da togliere per andare a regime
n=480;%input('# of measurements = ');%numero dati effettivi per la stima (50,300)
mordn=30;%input('model order = ');if isempty(mord);mord=10;end%ordine ARMAX (10,20)
off=1;%minimo ordine di modello
inp=input('1-WN 2-LP: ');if isempty(inp);inp=1;end
lowpar=0.8;%parametro per l'ingresso 2
al=200;%allunga con zeri le i.r.
nn=n+p+p2;
nt=500;%dimensione test set
frac=1/3% fraction of data used for validation/selection
ntval=floor(frac*(n+p));%dimensione validation set
red=n;
SV=5;%20
m=15;%15numero ingressi incogniti
mPEM=min(15,ceil((9/10)*(n+p)/(m+2)))%20;%max ordine PEM
Multi=[1:20];
nM=size(Multi,2);
ou=1;
oee=1;
BL=[];%i blocchi per la marginale, se [] disattivato

maxp_GLARS = 30 % maximum order for GLARS;

%m=4, n=170 o 270 o 370, p=30
%m=9, n=270, p=30
%m=14, n=370, p=30
%m=19, n=380, p=20, off=10, mordn=25


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% VARIABLES TO BE STORED %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
errnp0a=[];
errnp1a=[];
errnp0b=[];
errnp1b=[];
errnp0c=[];
errnp1c=[];
errbic0=[];
errbic1=[];
erraic0=[];
erraic1=[];
errOr0=[];
errOr1=[];
errgl0=[];
errgl1=[];
errgls0=[];
errgls1=[];


for jj=1:niter;jj
    clear numG numH den b U
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%GENERA MODELLO, TRAINING E TEST SET%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %T=1: modello normale
    %T=2: output error, no regressione
    %T=3: sparsita' casuale
    T=3;
    [den,numG,numH,mord,Ind(jj,:),V]=genmod(m,mordn,off,SV,T);
    Ind%disp(sqrt([V1 V2 V3 V4 V5]))
    [gtrue,htrue,a,b,y,U,Pr,ytest_tot,ytest,Utest,Ttest]=genri(den,numG,numH,inp,ou,oee,p,p2,nn,al,nt,lowpar,ntval);
    ya=y(1:(end-ntval+pRKHS),:);
    Ua=U(1:(end-ntval+pRKHS),:);
    yval=y((end-ntval+1):end,:);
    Uval=U((end-ntval+1):end,:);
    
    yaLARS=y(1:(end-ntval),:);
    UaLARS=U(1:(end-ntval),:);
    yvalLARS=y((end-ntval+1):end,:);
    UvalLARS=U((end-ntval+1):end,:);
    
    VarTest=sum((ytest-mean(ytest)).^2);
    
%     for i=1:m;
%         MOD(i,:,jj)=[numG(i,:) zeros(1,30-mord)];
%     end;
%     MOD(m+1,:,jj)=[den zeros(1,29-mord)];
%     MOD(m+2,:,jj)=[numH zeros(1,29-mord)];
%     DAT(:,:,jj)=[y U];%close all;plot(gtrue(1:100,1),'ko--');pause;close
%     all;plot(htrue(1:100,1),'ko--');pause
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %%%%%%%%%%%%%%%%%%%%MULTI-STEP OPTIMAL%%%%%%%%%%%%%%%%%%%%%%
%     RMSEbest=norm(Ttest*Pr-ytest)/sqrt(nt)
%     M = idpoly(1,gtrue(1:130,:)',htrue(1:130)',1,ones(m,1),1,1);
%     for gi=1:nM;
%         pred=predict(M,[ytest_tot Utest],Multi(gi));pred=iddata(pred);pred=pred.y;
%         RoptM(jj,gi)=norm(pred((p2+1):(p2+nt))-ytest)/sqrt(nt);
%         Ropt2M(jj,gi)=1-sum([pred((p2+1):(p2+nt))-ytest].^2)/VarTest;%norm([pred((p2+1):(p2+nt))-ytest]./ytest,1)/nt;
%     end
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%% STIMA NONPARAMETRICA %%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     lab='yy';%'yy' se c'e' parte autoregressiva e ingresso
%     mv=0;%mv = 1 se si stima una sola varianza, se no si stimano m varianze
%     mb=1;%mb = 1 se si stima un solo cambiamento di scala, se no si stimano m cambiamenti di scala
%     cn=2;%max numero poli negativi
%     ipn=-0.5*ones(1,cn);
%     LP=1;%se diverso da 1 non c'e il prior Laplace
%     [eta,Mnpa,Da,Sa,objval,A,iper,MP,nv]=SSpline(y,U,p,lab,mv,mb,ipn,[],LP,cn,red,BL);
%     %mv=1;LP=0;%una sola varianza ora, LP=1 sarebbe assurdo in questo caso
%     %[eta,Mnpb,Db,Sb,objval,A,iper,MP,nv]=SSpline(y,U,p,lab,mv,mb,ipn,[],LP,cn,red,BL);
%     [errnp0a,errnp1a,Errnpa(jj),RnpMa(jj,:),Rnp2Ma(jj,:),Da,Sa]=CompErr(Mnpa,gtrue,htrue,Ind(jj,:),errnp0a,errnp1a,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
%     %[errnp0b,errnp1b,Errnpb(jj),RnpMb(jj,:),Rnp2Mb(jj,:),Db,Sb]=CompErr(Mnpb,gtrue,htrue,Ind(jj,:),errnp0b,errnp1b,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
%     %HIP(jj,:)=ipera;
% 
%     clear AIK
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%% STIMA NONPARAMETRICA 2%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [Mnpb,beta,sigma,NT(jj),Mnpc,beta2,sigma2,NT2(jj)] = RKHS_Pred_LARS_Validation_all_k_v2(ya,Ua,1,30,yval,Uval,pRKHS);
    [errnp0b,errnp1b,Errnpb(jj),RnpMb(jj,:),Rnp2Mb(jj,:),Db,Sb]=CompErr(Mnpb,gtrue,htrue,Ind(jj,:),errnp0b,errnp1b,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
    [errnp0c,errnp1c,Errnpc(jj),RnpMc(jj,:),Rnp2Mc(jj,:),Dc,Sc]=CompErr(Mnpc,gtrue,htrue,Ind(jj,:),errnp0c,errnp1c,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
%     
% %      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %     %%%%%%%%%%%%%%%%%%%%% STIMA GLARS CP %%%%%%%%%%%%%%%%%%%%%%%%
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      dat_val=iddata(yval,Uval);
%     t=0;
%     for pl=1:1:maxp_GLARS;
%         t=t+1;
%         ll(t,:)=pl;
%         [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,pl,1,1,yvalLARS,[],UvalLARS);
%         predd=predict(Mgl,dat_val);predd=iddata(predd);predd=predd.y;
%         RM1(t)=sum((predd-yval).^2);
%         predd=predict(Mgls,dat_val);predd=iddata(predd);predd=predd.y;
%         RM2(t)=sum((predd-yval).^2);
%     end
%     [q,w1]=min(RM1);ll(w1,:)
%     [Mgl,Mgls] = Dyn_GLARS(ya,[],Ua,ll(w1,:),1);
%     [errnp0b,errnp1b,Errnpb(jj),RnpMb(jj,:),Rnp2Mb(jj,:),Db,Sb]=CompErr(Mgl,gtrue,htrue,Ind(jj,:),errgl0,errgl1,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);

    
    
    %data_test=iddata(ytest_tot,Utest);
    %yh=predict(Mgl,data_test);
    %eh=ytest_tot-yh.y;
    %err1=norm(eh)/sqrt(length(eh))
    %RnpMb_bis(jj,:)=norm(eh(pRKHS+1:end))/sqrt(length(eh(pRKHS+1:end)))
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%% STIMA GLARS VALIDATION %%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    dat_val=iddata(yval,Uval);
    t=0;
    for pl=1:1:maxp_GLARS;
        t=t+1;
        ll(t,:)=pl;
        [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,pl,1,0,yvalLARS,[],UvalLARS);
        predd=predict(Mgl,dat_val);predd=iddata(predd);predd=predd.y;
        RM1(t)=sum((predd-yval).^2);
        predd=predict(Mgls,dat_val);predd=iddata(predd);predd=predd.y;
        RM2(t)=sum((predd-yval).^2);
    end
    [q,w1]=min(RM1);ll(w1,:)
    [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,ll(w1,:),1,0,yvalLARS,[],UvalLARS);
    [errgl0,errgl1,Errgl(jj),RglM(jj,:),Rgl2M(jj,:),Da,Sa]=CompErr(Mgl,gtrue,htrue,Ind(jj,:),errgl0,errgl1,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
    
%     [q,w2]=min(RM2);ll(w2,:)
%     [Mgl,Mgls] = Dyn_GLARS(ya,[],Ua,ll(w2,:),1);
%     [errgls0,errgls1,Errgls(jj),RglsM(jj,:),Rgls2M(jj,:),Da,Sa]=CompErr(Mgls,gtrue,htrue,Ind(jj,:),errgls0,errgls1,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%STIMA PEM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    dat=iddata(ya,Ua);
    dat_test=iddata(ytest_tot,Utest);
    t=0;
    for i1=1:mPEM;
            i2=i1*ones(1,m);
            i3=i1;
            t=t+1;
            %tic
            Mpem=pem(dat,[i1 i2 i3 zeros(1,m+1) ones(1,m)],'LimitError',0);
            %toc
            %pause
            %tic
            %Mpem=pem(dat,i1);
            %toc 
            %pause
            predd=predict(Mpem,dat_val);predd=iddata(predd);predd=predd.y;
            RM(t)=sum((predd-yval).^2);
            ll1(t,:)=[i1 i2 i3];
    end
    [q,w1]=min(RM);ll1(w1,:)
    Mpem=pem(dat,[ll1(w1,:) zeros(1,m+1) ones(1,m)],'LimitError',0);
    [errbic0,errbic1,Errbic(jj),RbicM(jj,:),Rbic2M(jj,:),Dbic,Sbic]=CompErr(Mpem,gtrue,htrue,Ind(jj,:),errbic0,errbic1,ytest_tot,ytest,Utest,p2,nt,VarTest,Multi);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%RESULTS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %close all
    figure(1)
    hold off
    subplot(2,1,1);
    plot(mean(Rbic2M),'b*');hold on;
    %plot(mean(Rnp2Ma),'ko');
    plot(mean(Rgl2M),'go');
    %plot(mean(Rgls2M),'g*');
    plot(mean(Rnp2Mb),'kd');
%    plot(mean(Rnp2Mc),'ko');
    %plot(mean(ROr2M),'r.');
    %plot(mean(Ropt2M),'r.');

    subplot(2,1,2);
    hold off
    plot(mean(RbicM),'b*');hold on
    plot(mean(RglM),'go');
    %plot(mean(RglsM),'g*');
    plot(mean(RnpMb),'kd');
 %   plot(mean(RnpMc),'ko');
    %plot(mean(RnpMb),'go');
    %plot(mean(ROrM),'r.');
    %plot(mean(RoptM),'r.');
    pause(1)
    
    %disp('np')
    %[mean(errnp0a) mean(errnp1a) mean(Errnpa)  mean(mean(Rnp2Ma))]%norma di quelle che dovrebbero essere a zero, errore relativo sulle norme delle altre, 
                                           %errore relativo sulle norme
                                           %generale (a differenza degli altri, un solo indice per corsa)
    
                                
                                     
    disp('GLARS CP')                                       
    [mean(errnp0b) mean(errnp1b) mean(Errnpb)  mean(mean(Rnp2Mb))  mean(RnpMb(:,1)) ]
    
   
%     disp('np3')                                       
%     [mean(errnp0c) mean(errnp1c) mean(Errnpc)  mean(mean(Rnp2Mc))  mean(RnpMc(:,1))]
    
    disp('GLARS val')
    [mean(errgl0) mean(errgl1) mean(Errgl)  mean(mean(Rgl2M))  mean(RglM(:,1))]
    
%     disp('GLARSls')
%     [mean(errgls0) mean(errgls1) mean(Errgls)  mean(mean(Rgls2M))  mean(RglsM(:,1))]
%     
    disp('PEM+Val')
    [mean(errbic0) mean(errbic1) mean(Errbic)  mean(mean(Rbic2M))  mean(RbicM(:,1))]
    
    %disp('PEM+AICC')
    %[mean(erraic0) mean(erraic1) mean(Erraic)  mean(mean(Raic2M))]
    
    %disp('PEM+OR')
    %[mean(errOr0) mean(errOr1) mean(ErrOr)  mean(mean(ROr2M))]
    
    %boxplot([errnp0' errbic0'])%sparsita'
    %bode(M,'r',Mnp,'k',Mpem,'b')
lab={'GLcp' 'GLVAL'  'PEM'};
%close all
w=max(size(RbicM,1));
if w>1
    figure(2)
    hold off %boxplot([ RnpMb(1:w,1)  RnpMc(1:w,1)  RglM(1:w,1)  RglsM(1:w,1)  RbicM(1:w,1)],'labels',lab)
    boxplot([ RnpMb(1:w,1)   RglM(1:w,1)    RbicM(1:w,1)],'labels',lab)
    drawnow
%     if w==2 
%         pause
%     end;
    pause(5)
    
%     if w==2
%         pause
%     end;

end;

end


% lab={'Np1' 'Np2' 'Np3' 'GL1' 'GL2' 'PEM'};
% close all
% w=max(size(RbicM));
% figure(1)
% boxplot([RnpMa(1:w,1)  RnpMb(1:w,1)  RnpMc(1:w,1)  RglM(1:w,1)  RglsM(1:w,1)  RbicM(1:w,1)],'labels',lab)
% w=max(size(errbic1));
% figure(2)
% boxplot([errnp1a(1:w)'  errnp1b(1:w)'   errnp1c(1:w)'   errgl1(1:w)'  errgls1(1:w)' errbic1(1:w)'],'labels',lab)
% [sum(errnp0a==0)/max(size(errnp0a))  sum(errnp0b==0)/max(size(errnp0b))  sum(errnp0c==0)/max(size(errnp0c))  sum(errgl0==0)/max(size(errgl0))  sum(errgls0==0)/max(size(errgls0))...
%     sum(errbic0==0)/max(size(errbic0))]


lab={'GLcp'  'GLval'  'PEM'};
close all
w=max(size(RbicM));
figure(1)
%boxplot([ RnpMb(1:w,1)  RnpMc(1:w,1)  RglM(1:w,1)  RglsM(1:w,1)  RbicM(1:w,1)],'labels',lab)
boxplot([ RnpMb(1:w,1)   RglM(1:w,1)    RbicM(1:w,1)],'labels',lab)
w=max(size(errbic1));
figure(2)
%boxplot([ errnp1b(1:w)'   errnp1c(1:w)'   errgl1(1:w)'  errgls1(1:w)' errbic1(1:w)'],'labels',lab)
%boxplot([ errnp1b(1:w)'   errnp1c(1:w)'   errgl1(1:w)'   errbic1(1:w)'],'labels',lab)
boxplot([Errnpb'     Errgl'   Errbic'],'labels',lab)
%[sum(errnp0b==0)/max(size(errnp0b))  sum(errnp0c==0)/max(size(errnp0c))  sum(errgl0==0)/max(size(errgl0))  sum(errgls0==0)/max(size(errgls0))...
%    sum(errbic0==0)/max(size(errbic0))]

[sum(errnp0b==0)/max(size(errnp0b))    sum(errgl0==0)/max(size(errgl0)) ...
    sum(errbic0==0)/max(size(errbic0))]


% errnp1a
% errbic1
% close all
% plot([1:100],htrue(1:100),'r',[1:100],Sa(1:100),'k',[1:100],Sbic(1:100),'b')
% pause
% for i=1:size(gtrue,2);
%     plot([1:100],gtrue(1:100,i),'r',[1:100],Da(1:100,i),'k',[1:100],Dbic(
%     1:100,i),'b')
%     pause
% end


%%%%%%%%%%%%%%%%%%%%%%

% ONLY FOR DISPLYING FINAL RESULTS

disp('np2')                                       
    [mean(errnp0b) mean(errnp1b) mean(Errnpb)  mean(mean(Rnp2Mb))  mean(RnpMb(:,1)) mean(RnpMb_bis(:,1))]
    
   
    % disp('np3')                                       
    % [mean(errnp0c mean(errnp1c) mean(Errnpc)  mean(mean(Rnp2Mc))  mean(RnpMc(:,1))]
    
    disp('GLARS')
    [mean(errgl0) mean(errgl1) mean(Errgl)  mean(mean(Rgl2M))  mean(RglM(:,1))]
    
%     disp('GLARSls')
%     [mean(errgls0) mean(errgls1) mean(Errgls)  mean(mean(Rgls2M))  mean(RglsM(:,1))]
%     
    disp('PEM+Val')
    [mean(errbic0) mean(errbic1) mean(Errbic)  mean(mean(Rbic2M))  mean(RbicM(:,1))]
    
    
    figure(3) 
plot(mean(Rnp2Mb),'o')
hold on
plot(mean(Rgl2M),'+')
plot(mean(Rbic2M),'*')
title('COD, k-step-ahead prediction')
xlabel('prediction horizon k')
Legend('SS-GLAR','VARX-GLAR','PEM')

figure(4) 
plot(mean(RnpMb),'o')
hold on
plot(mean(RglM),'+')
plot(mean(RbicM),'*')
title('k-step-ahead RMS error')
xlabel('prediction horizon k')
Legend('SS-GLAR','VARX-GLAR','PEM')

figure(5)
subplot(211),plot(mean(Rnp2Mb),'o')
hold on
plot(mean(Rgl2M),'+')
plot(mean(Rbic2M),'*')
title('average COD, k-step-ahead prediction')
xlabel('prediction horizon k')
Legend('SS-GLAR','VARX-GLAR','PEM')

subplot(212),plot(mean(RnpMb),'o')
hold on
plot(mean(RglM),'+')
plot(mean(RbicM),'*')
title('average k-step-ahead RMS error')
xlabel('prediction horizon k')
Legend('SS-GLAR','VARX-GLAR','PEM')




lab={'SS-GLAR' 'VARX-GLAR'  'PEM'};

w=max(size(RbicM));
figure(1)
%boxplot([ RnpMb(1:w,1)  RnpMc(1:w,1)  RglM(1:w,1)  RglsM(1:w,1)  RbicM(1:w,1)],'labels',lab)
boxplot([ RnpMb(1:w,1)  RglM(1:w,1)      RbicM(1:w,1)],'labels',lab)
w=max(size(errbic1));
figure(2)
%boxplot([ errnp1b(1:w)'   errnp1c(1:w)'   errgl1(1:w)'  errgls1(1:w)' errbic1(1:w)'],'labels',lab)
%boxplot([ errnp1b(1:w)'   errnp1c(1:w)'   errgl1(1:w)'   errbic1(1:w)'],'labels',lab)
boxplot([Errnpb'   Errnpc'   Errgl'   Errbic'],'labels',lab)

