%
%  fit radial basis function net to one dimensional
%  input-output data
%
%   [ yhat, cc ] = fitrbfn( u, y, k )
%
%  yhat  = approximated output data
%  cc    = net parameters = [ c tk sigma ]
%  c     = rbf coefficients
%  tk    = rbf centres
%  sigma = rbf widths
%  u     = input data
%  y     = output data
%  k     = number of nodes (rbf's)
%
function [ yhat, cc ] = fitrbfn( u, y, k )
%
% use clustering to locate basis functions
[ tk, alloc ] = kmeans( u(:), k );
%
% use nearest neighbour to size basis functions
for i = 1:k,
  d = [];
  for j = 1:k,
    if j ~= i,
      d = [ d sqrt(abs(tk(i)*tk(i)-tk(j)*tk(j))) ];
    end;
  end;
  d = sort( d );
  sigma(i) = sqrt(d(1)*d(2));
end;
%
% construct regression matrix
n = length( u );
for i = 1:n,
  for j = 1:k,
    a(i,j) = exp( -((u(i)-tk(j))./sigma(j)).^2 );
  end;
end;
b = y(:); 
%
% calculate coefficients
c = a\b;
%
% calculate approximation
yhat = a * c;
sd = sqrt( ( y - yhat )' * ( y - yhat ) / n );
sd
%
cc = [ c tk sigma(:) ];
%

