clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imposta i parametri del modello generatore
% dei dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('seed',54321);
rand('seed',54321);
N=400;                     % Numerosità campionaria
ntrain = floor(2*N/3);     % Numerosità del training set
ninp = 5;                  % No. variabili X di input
devstd = 0.2;              % Deviazione standard della componente stocastica
goal = 0;                  % Errore medio quadratico da raggiungere
spread = 1;                % Parametro di scala sigma
Nhid_Range = [1:4 5:5:80]; % Gamma di valori per il numero massimo di neuroni

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = devstd*randn(N,1);                      % Componente stocastica
y = ydet + eps;                               % Variabile dipendente y

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Costruisci training and testing set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dati non riscalati
X_train    = X(1:ntrain,:);
y_train    = y(1:ntrain);
ydet_train = ydet(1:ntrain);
X_test     = X(ntrain+1:end,:);
y_test     = y(ntrain+1:end);
ydet_test  = ydet(ntrain+1:end);
% Dati riscalati
Xtrain = Xn(:,1:ntrain);
ytrain = yn(1:ntrain);
Xtest  = Xn(:,ntrain+1:end);
ytest  = yn(ntrain+1:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo dell'errore in-sample ed out-of-sample
% per diversi valori del numero massimo di neuroni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UT_train = zeros(length(Nhid_Range),1);
UT_test = UT_train;
no_iter = length(Nhid_Range);
for k = 1:no_iter
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Costruzione e stima rete RBF
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    MaxNeu = Nhid_Range(k);
    optnet=newrb(Xtrain,ytrain,goal,spread,MaxNeu,500);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Previsione con la rete neurale
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    yout_train=sim(optnet, Xtrain);
    yout_test=sim(optnet, Xtest);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Postprocessamento output
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if verLessThan('matlab', '7.1.0')
        yhat_train = postmnmx(yout_train,miny,maxy);
        yhat_test  = postmnmx(yout_test ,miny,maxy);
    else
        yhat_train = mapminmax('reverse',yout_train,psy);
        yhat_test  = mapminmax('reverse',yout_test ,psy);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calcolo della bontà di adattamento
    % ai dati mediante indice U-Theil
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    UT_train(k) = 100*sum((y_train-yhat_train').^2)/sum((y_train-mean(y_train)).^2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calcolo della capacità predittiva della rete
    % mediante indice U-Theil
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    UT_test(k) = 100*sum((y_test-yhat_test').^2)/sum((y_test-mean(y_test)).^2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rappresentazione grafica della bontà della 
% rete di adattarsi ai dati (goodness of fitting) 
% e della sua capacità previsiva
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(Nhid_Range,UT_train,'k-',Nhid_Range,UT_test,'k:');
xlabel('Numero neuroni strato nascosto')
ylabel('U-Theil')
legend('In-sample','Out-of-sample');

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