introexample.m

人工神经网络的源码编程

%%
%% Example for introductory (tutorial) manual pages.
%% Fits a 1-dimensional sine wave.
%%

%%
%% Prepare figures.
%%

figure(1)
pos = get(1, 'Position');
set(1, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample1', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(2)
pos = get(2, 'Position');
set(2, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample3', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(3)
pos = get(3, 'Position');
set(3, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample4', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(4)
pos = get(4, 'Position');
set(4, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample5', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(5)
pos = get(5, 'Position');
set(5, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample6', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(6)
pos = get(6, 'Position');
set(6, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample7', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(7)
pos = get(7, 'Position');
set(7, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample8', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])
figure(8)
pos = get(8, 'Position');
set(8, ...
  'Position', [pos(1) pos(2) 400 400], ...
  'NumberTitle', 'off', ...
  'Name', 'introExample9', ...
  'PaperType', 'a4letter', ...
  'InvertHardCopy', 'on', ...
  'PaperPosition', [0.5 0.5 4 4])

%%
%% Get training data and plot it.
%%

p = 50;
sigma = 0.2;
NEW = 0;
if NEW
  x = 1 - rand(1,p);
  y = sin(10 * x)' + sigma * randn(p,1);
  save introExample x y
else
  load introExample
end
figure(1)
hold off
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
title('Figure 1')
drawnow

%%
%% Get test data and plot target.
%%

pt = 500;
xt = linspace(0,1,pt);
yt = sin(10 * xt)';
hold on
plot(xt, yt, 'b--')

%%
%% Plot single Gaussian radial basis function.
%%

cj = 0.5;
rj = 0.1;
ht = exp(-(xt - cj).^2 / rj^2);
figure(2)
hold off
plot(xt, ht, 'g-')
axis([0 1 -0.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [0 1])
title('Figure 3')
drawnow

%%
%% Design matrices.
%%

c = x;
m = p;
r = 0.1;
H = rbfDesign(x, c, r, 'c');
Ht = rbfDesign(xt, c, r, 'c');

%%
%% Train with normal equation.
%%

w = inv(H' * H) * H' * y;

%%
%% Plot the output.
%%

ft = Ht * w;
figure(3)
hold off
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
hold on
plot(xt, yt, 'b--')
plot(xt, ft, 'g-')
title('Figure 4')
drawnow

%%
%% Forward selection: predict and plot.
%%

subset = forwardSelect(H, y, '-t BIC -v');
Hs = H(:,subset);
w = inv(Hs' * Hs) * Hs' * y;
ft = Ht(:,subset) * w;
figure(4)
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
hold on
plot(xt, yt, 'b--')
plot(xt, ft, 'g-')
title('Figure 5')
drawnow

%%
%% Mean square error.
%%

(ft - yt)' * (ft - yt) / pt

%%
%% Global ridge with guessed regularisation paremeter.
%%

lambda = 1e-4;
w = inv(H' * H + lambda * eye(m)) * H' * y;
ft = Ht * w;
figure(5)
hold off
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
hold on
plot(xt, yt, 'b--')
plot(xt, ft, 'g-')
title('Figure 6')
drawnow

%%
%% Mean square error.
%%

(ft - yt)' * (ft - yt) / pt

%%
%% Actively choose the right regularisation parameter.
%%

lambda = globalRidge(H, y, 0.1, 'BIC -v')

%%
%% Global ridge with optimal lambda.
%%

w = inv(H' * H + lambda * eye(m)) * H' * y;
ft = Ht * w;
figure(6)
hold off
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
hold on
plot(xt, yt, 'b--')
plot(xt, ft, 'g-')
title('Figure 7')
drawnow

%%
%% Mean square error.
%%

(ft - yt)' * (ft - yt) / pt

%%
%% Explicit calculation and plot of BIC for range of lambda.
%%

lambdas = logspace(-4, 2, 100);
bics = zeros(1,100);
b = 0;
for i = 1:100
  b = overWrite(b, i);
  bics(i) = predictError(H, y, lambdas(i), 'BIC');
end
overWrite(b);
bic = predictError(H, y, lambda, 'BIC');
figure(7)
hold off
loglog(lambdas, bics, 'm-')
hold on
loglog(lambda, bic, 'r*')
title('Figure 8')
drawnow

%%
%% Try local ridge regression.
%%

lambdas = localRidge(H, y, lambda, '-v');
finite = find(lambdas ~= Inf)

%%
%% Work out fit and plot it.
%%

Hl = H(:,finite);
lambdas = lambdas(finite);
w = inv(Hl' * Hl + diag(lambdas)) * Hl' * y;
ft = Ht(:,finite) * w;
figure(8)
hold off
plot(x, y, 'c+')
axis([0 1 -1.5 1.5])
set(gca, 'XTick', [0 0.5 1])
set(gca, 'YTick', [-1 0 1])
hold on
plot(xt, yt, 'b--')
plot(xt, ft, 'g-')
title('Figure 9')
drawnow

%%
%% Mean square error.
%%

(ft - yt)' * (ft - yt) / pt