% 1. Create the 6 trn/val/tst input/target subsets
% 2. For a c-class classifier, the columns of the target matrices
%    are all columns of the c-dimensional unit matrix eye(c) with
%    the single 1 indicating the corresponding class. Inputs will
%    be assigned to the class coresponding to the largest output.
% 3. Standardize (e.g., mapstd) the trn input and normalize the val
and
%    trn inputs with the means and variances from trn.
% 4. It is sufficient to use a single hidden layer so that the node
%    topology for I-dimensional inputs and O-dimensional outputs
%    is I-H-O (O=c);
% 5. Calculate the trn/val/tst MSEs for the naive constant models.
%    For example,
%    ytrn00          = repmat(mean(ttrn,2),1,Ntrn);
%    MSEtrn00    = mse(ttrn-ytrn00)
%    Neq00         = Ntrn*O                  % No. of equations
%    Nw00          = O                          % No. of weights
%    MSEtrn00a = sse(ttrn-ytrn00)/(Neq00-Nw00)  % MSE "a"djusted to
mitigate the optimistic bias caused by estimating generalization
performance using the same data that was used for weight estimation.
%    etc
% 6. The number of training equations is Neq = Ntrn*O and the number
%    of unknown weights is Nw = (I+1)*H+(H+1)*O. For training to
%    convergence, a necessary criterion for a unique solution is
%    Neq >= Nw or H <= Hub where the upperbound is given by
%    Hub = (Neq-O)/(I+O+1). However, weight estimates are more
%    resistant to noise and measurement errors when H << Hub. The
%    optimal value, Hopt, is usually estimated via trial and error.
% 7. To avoid nonoptimal local min solutions, choose the best of
%    many designs obtained by using multiple random weight
%    intializations for each candidate value of H.
% 8. Determine a search grid for H: (e.g., Hmin:dH:Hmax) and the
%    number of random weight initialization trials, Ntrials, for
%    each H candidate.
% 9. Initialize the RAND generator and create numH*Ntrials candidate
%    designs using an outer loop of length numH = length(Hmin:dH:Hmax)
%    and an inner loop of length Ntrials.
% 10. The canonical Bayesian classification goal is to minimize the
%     total misclassification risk, given class a priori
probabilities
%     and misclassification costs. However, most practical approaches
%     1. Use duplication or simulation of additional samples to
prevent
%        significantly unbalanced class mixtures of training data.
%     2. Minimize the mean-square error of the difference between the
%        output estimates of posterior probability and the target
%        matrices.
%     3. Use the estimates of val set mixture misclassification rate
%        to select the best group of designs.
%     4. Either choose the best design or combine the outputs from
%        a committee of best designs.
%
%   Consider the following demo of an I-H-O Elliptical Basis
Function
%   (EBF) design using NEWFF and RADBAS. The search is over H =
2.^(0:6),
%   Ntrials = 10. The training goal is an adjusted training R^2 of
0.99
%   which is obtained half the time when H = 32 .
%   Modification of the code to include a validation set to choose
the
%   best designs and a test set to estimate performance on unseen
data
%   is left to the reader.


  close all, clear all, clc, k=0

  N1 = 30, N2 = 30, N3 = 30, N4 = 30
  randn('state',4151941)
  p1 = repmat([-1;-1],1,N1)+0.6*randn(2,N1);
  p2 = repmat([ 1;-1],1,N2)+0.6*randn(2,N2);
  p3 = repmat([ 1; 1],1,N3)+0.6*randn(2,N3);
  p4 = repmat([-1; 1],1,N4)+0.6*randn(2,N4);

  eye4 = eye(4);
  t1 = repmat(eye4(:,1),1,N1)
  t2 = repmat(eye4(:,2),1,N2)
  t3 = repmat(eye4(:,3),1,N3)
  t4 = repmat(eye4(:,4),1,N4)

  k=k+1,figure(k),hold on
  plot(-3:3,zeros(1,7),'k')
  plot(zeros(1,7),-3:3,'k')
  plot(p1(1,:),p1(2,:),'ko')
  plot(p2(1,:),p2(2,:),'ro')
  plot(p3(1,:),p3(2,:),'bo')
  plot(p4(1,:),p4(2,:),'go')

  ptrn = [p1 p2 p3 p4];
  ttrn = [t1 t2 t3 t4];
  [I Ntrn] = size(ptrn)         % [2 120]
  [O Ntrn] = size(ttrn)         % [4 120]
  Neq      = Ntrn*O             % 480

  MSEtrn00  =  mse(ttrn - repmat(mean(ttrn,2),1,Ntrn)) % 0.1875
  MSEtrn00a =  mean(var(ttrn'))                        % 0.1891

%   Nw00      =  O % No. of weights
%   MSEtrn00  =  sse(ttrn - repmat(mean(ttrn,2),1,Ntrn))/Neq        %
0.1875
%   MSEtrn00a =  sse(ttrn - repmat(mean(ttrn,2),1,Ntrn))/(Neq-Nw00) %
0.1891

  Hub  = (Neq-O)/(I+O+1)    % 68

  Hmin = 16
  dH   = 4   % Minimum allowable value is 1
  Hmax = 34

  Ntrials = 10
  rand('state',0)

  for i = 1:7   %H = [1 2 4 8 16 32 64]
      H = 2^(i-1)
      Nw = (I+1)*H+(H+1)*O
      for j = 1:Ntrials
          net = newff(minmax(ptrn),[H O],{'radbas', 'purelin'});
          net.trainParam.goal   = 0.01*(Neq-Nw)*MSEtrn00/Neq; %R2trna
= 0.99
          net.trainParam.show   = inf;
%           net.trainParam.epochs = 500;
          [net tr] = train(net,ptrn,ttrn);
          Nepochs(j,i) = tr.epoch(end);
          MSEtrn       = tr.perf(end);
          R2trn(j,i)   = 1 - MSEtrn/MSEtrn00;
          R2trna(j,i)  = 1- (Neq/(Neq-Nw))*MSEtrn/MSEtrn00;
%           MSEval       = mse(tval-sim(net,pval));
%           R2val(j,i)   = 1- MSEval/MSEval00;
%           MSEtst       = mse(ttst-sim(net,ptst));
%           R2tst(j,i)   = 1- MSEtst/MSEtst00;
      end
  end

   H = 2.^(0:6)
   Nepochs
   R2trn
   R2trna

      H =    1     2     4     8    16    32    64

Nepochs =   60    78    53   500   500   500    58
            62    37    41   468   359   500    50
            76    49   500   209   500    82    46
            91    31    90   155   500   454    56
            58    32    82   500   500   198    54
            61    31    49   500   500   500    44
            30    54    58   145   500   500    52
            70    53   500   141   500   237    40
            58    56    49   500   500   115    62
            61    49   143   239   500   500    56

R2trn = 0.1780    0.5476    0.8271    0.8685    0.9362    0.9932
0.9994
        0.2207    0.5696    0.7972    0.8737    0.9304    0.9867
0.9994
        0.2930    0.5669    0.7998    0.8784    0.9517    0.9948
0.9994
        0.2781    0.5673    0.8131    0.8631    0.9295    0.9948
0.9994
        0.1780    0.5669    0.8124    0.8775    0.9383    0.9950
0.9994
        0.2752    0.5673    0.8271    0.8576    0.9421    0.9937
0.9994
        0.2915    0.5669    0.8310    0.8924    0.9324    0.9853
0.9994
        0.1780    0.5649    0.8129    0.8719    0.9291    0.9948
0.9995
        0.2207    0.4934    0.7899    0.8817    0.9659    0.9948
0.9994
        0.2207    0.5673    0.7992    0.8682    0.9269    0.9713
0.9994

R2trna = 0.1588    0.5300    0.8148    0.8498    0.9158    0.9871
0.9903
         0.2024    0.5529    0.7827    0.8557    0.9083    0.9746
0.9902
         0.2764    0.5501    0.7855    0.8610    0.9363   *0.9901
0.9900
         0.2612    0.5505    0.7998    0.8435    0.9070   *0.9900
0.9904
         0.1588    0.5501    0.7990    0.8601    0.9186   *0.9904
0.9901
         0.2582    0.5505    0.8148    0.8373    0.9237    0.9881
0.9904
         0.2749    0.5501    0.8189    0.8770    0.9109    0.9721
0.9903
         0.1588    0.5480    0.7995    0.8536    0.9065   *0.9900
0.9908
         0.2024    0.4737    0.7749    0.8649    0.9551   *0.9901
0.9903
         0.2024    0.5505    0.7849    0.8494    0.9036    0.9454
0.9903