Elliptical Basis Function NN Design

[Greg Heath]

% 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