close all; clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imposta i parametri della simulazione
% e della rete neurale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('seed',16220326)
rand('seed',16220326)
N=200;    % Numerosità campionaria
ninp=5;   % No. variabili X di input
nhid=20;  % No. di neuroni nello strato nascosto
sigma=0.2 % Deviazione standard della componente stocastica

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Genera i dati campionari
% Data generating process: 
% y = 1.2*(x1^3 - 0.5) + sin(3.46*pi*x1) + eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X =rand(N,ninp);                              % Variabili indipendenti
ydet= 1.2*(X(:,1).^3 -.5)+sin(X(:,1)*pi*1.1); % Componente deterministica
eps = sigma*randn(N,1);                       % Componente stocastica
y = ydet + eps;                               % Variabile dipendente y

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rappresentazione grafica del problema di stima
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
hold on;
x1 = X(:,1);
[x1sort,order]=sort(x1);
plot(x1,y,'b.');                              % Variabile dipendente y
plot(x1sort,ydet(order),'r-','LineWidth',2);  % Componente deterministica
box on;
legend('Dati osservati','Comp.determinist.',...
       'location','SouthEast');
hold off;
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preprocessamento dei dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verLessThan('matlab', '7.1.0')
    [Xn,minX,maxX,yn,miny,maxy] = premnmx(X',y');
else
    [Xn,psx] = mapminmax(X',-1,1);
    [yn,psy] = mapminmax(y',-1,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Crea rete neurale feedforward
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
input_minmax = minmax(Xn);
net=newff(input_minmax,[nhid 1],{'tansig','tansig'},'trainlm');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imposta pesi della rete neurale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net=init(net);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Training della rete neurale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parametri per lo stop dell'algoritmo di ottimizzazione
% No. massimo di iterazioni dell'algoritmo di ottimizzazione
net.trainParam.epochs=500;
% Ferma l'apprendimento se l'errore sul training set
% scende al di sotto di questa soglia (0.001)
net.trainParam.goal = 10^(-3);
% Training della rete
optnet=train(net, Xn , yn); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Previsione con la rete neurale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yout=sim(optnet, Xn); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Postprocessamento output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verLessThan('matlab', '7.1.0')
    yhat = postmnmx(yout,miny,maxy);
else
    yhat = mapminmax('reverse',yout,psy);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo di un indicatore della bontà di 
% adattamento del modello ai dati: 
% l'indice U di Theil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yhat = yhat';
UT = 100*sum((y-yhat).^2)/sum((y-mean(y)).^2);
display(['U-Theil = ', num2str(UT), '%']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rappresentazione grafica della 
% capacita' predittiva della rete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
hold on;
plot(x1,y,'b.')
plot(x1sort,yhat(order),'-','Color',[0.23 0.62 0.17]);
plot(x1sort,ydet(order),'r-','LineWidth',2);
box on;
legend('Dati osservati','Previsioni rete','Comp.determinist.',...
       'location','SouthEast');
title(['U-Theil = ', num2str(UT), '%'])
hold off;
print regr_RNA_1.wmf -dmeta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analisi dei residui
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)
res1 = y - yhat;
% Grafico valori veri contro valori stimati dalla rete neurale
subplot(2,2,1)
plot(x1,y,'.',x1sort,yhat(order),'-')
title('Y vs. X e previsione modello')
title(['Y e previsione vs. X'])
% Grafico dei residui del modello contro il valore di x
subplot(2,2,2)
plot(x1sort,res1(order))
axis tight
xlabel('X')
ylabel('Residui')
title('Residui modello vs. valori di X')
% Istogramma dei residui e curva gaussiana per
% il confronto delle distribuzioni
subplot(2,2,3)
me1=mean(res1);
sd1=std(res1);
res1std=(res1-me1)/sd1;
[fr1,xc1]=hist(res1std,40);
delta1=xc1(2)-xc1(1);
fr1 = fr1/sum(fr1.*delta1);
% Trasformo le frequenze in densita' e traccio l'istogramma
bar(xc1,fr1,1);
hold on;
% Traccio il grafico della guassiana
xpdf=[-5:.05:5];
plot(xpdf,normpdf(xpdf),'r');
% Test di gaussianita' dei residui con il metodo di Jarque-Bera
warning off;
[h,p]=jbtest(res1);
title(['p-value JB test = ',num2str(round(1000*p)/1000)])
% Il Q-Q Plot dei residui del modello per valutarne la gaussianita'
subplot(2,2,4)
qqplot(res1)
print regr_RNA_2.wmf -dmeta

% Marco Sandri
% Università di Brescia
% Dipartimento Metodi Quantitativi
% info@msandri.it