batch = 1;  % for generating plots
if batch, rand('seed',0);  randn('seed',0);  end

xlim = [0 10];
xpv = xlim(1):0.01:xlim(2);
ypv = -0.1*xpv.^3 + 2*(xpv-3).^2 + 20;
np = 40;
xv = (xlim(2)-xlim(1))*rand(1,np) + xlim(1);
yv = -0.1*xv.^3 + 2*(xv-3).^2 + 20 + 1*randn(size(xv));

% Problem statement
fh = figure;
plot(xpv,ypv,'k-');
hold on;  sh = scatter(xv,yv,'r+');  hold off;
axis([0 10 0 40]);
title('Points from general curve');
exportfig(gcf,'rbf_example_fig01.eps','width',10,'height',6.5,'fontmode','fixed', 'fontsize',8,'Color','cmyk');

% Radial basis functions
nrbf = 8;
rbfp = (0:nrbf-1)*(xlim(2)-xlim(1))/(nrbf-1) + xlim(1);
rbfsig = 1;
rbfm = zeros(nrbf,length(xpv));
for j=1:nrbf
  rbfm(j,:) = exp(-1/(2*rbfsig)*(xpv-rbfp(j)).^2);
end
fh = figure;
plot(xpv,rbfm')
title('Radial basis functions');
exportfig(gcf,'rbf_example_fig02.eps','width',10,'height',6.5,'fontmode','fixed', 'fontsize',8,'Color','cmyk');

% Least squares solution
A = zeros(np,nrbf);
for j=1:nrbf
  A(:,j) = exp(-1/(2*rbfsig)*(xv'-rbfp(j)).^2);
end
theta = inv(A'*A)*A'*yv';

% Original and estimated curve
fh = figure;
rypv = rbfm'*theta;
crbfm = repmat(theta,1,size(rbfm,2)).*rbfm;
plot(xpv,ypv,'k-',xpv,rypv,'r-');
hold on;  sh = scatter(xv,yv,'r+');  hold off;
hold on;  ph = plot(xpv,crbfm');  hold off;
set(ph,'color','g','LineStyle','--')
axis([0 10 0 40]);
title('RBF regression curve');
exportfig(gcf,'rbf_example_fig03.eps','width',10,'height',6.5,'fontmode','fixed', 'fontsize',8,'Color','cmyk');
