Re: Automatic polynomial approximation

[heberf]

This is not exactly what you want but here is a function which fits a curve to 
data which is non-parametric, meaning that no functional form is imposed.  The 
procedure is called kernel regression and the function is closely related to 
kernel.m which I submitted to the sources list several years ago.  Kernel 
regression is simpler than splines and better than polynomial fitting.  In fact 
it can be shown that kernel regression is equivalent to locally fitting a 
quadratic function to your data.

A good book on non-parametric regression in general and kernels in particular is

@Book{Hardle,
  author =       {Hardle, Wolfgang},
  title =        {Applied nonparametric regression},
  publisher =    {Cambridge University Press},
  year =         {1990},
  OPTaddress =   {New York},
}



function [xx, yy] = kreg(x,y,n,c)

usage [xx, yy] = kreg(x,y,n,c)

# Curve fitting via kernel regression. This version uses a Gaussian kernel
# and a bandwidth equal to 
#
#      h(m) = c*std(x)*m^(.2)
#
# where std(x) is the standard deviation of the data and m is the number of
# observations in x. 
#
# Given a third argument, use that as the number of evenly spaced points
# where the density will be evaluated.  Default is 40.
# The fourth argument,c, is a smoothing parameter.  Larger values of c lead to 
#  smoother curves.  The default value of c is one.
#
# With no output arguments the curve is plotted.  Given two output
# arguments the values xx and yy are returned such that plot(xx,yy) will
# plot the curve.  

# author: Heber Farnsworth 
# Dept of Finance,
# Olin School of Business
# Washington University in St. Louis
# address@hidden  


  if (nargin < 1 || nargin > 4) 
    usage ("[xx, yy] = kreg (x, y, n, c)");
  endif

  if nargin < 4
    c = 1;
  elseif c <=0 
    error ("kreg: second argument must be positive");
  endif
 
  if nargin < 3 
    n = 40;
  elseif n <= 2
    error ("kreg: number of points must be greater than 2");
  endif
  
 if nargin < 2 
   error ("kreg: x and y data vectors required");
 elseif (is_vector (x) && is_vector (y))
    max_val = max (x);
    min_val = min (x);
    m = max(size(x));
    my = max(size(y));
    if m != my 
      error ("kreg: data vectors must have same length")
    endif
    if m > size(x,1)
      x = x';
    endif
    if my > size(y,1)
      y = y';
    endif
  else
    error ("kreg: first two arguments must be vectors");
  endif

  s = std(x);
  h = c*s/(m^(.2));
  z = linspace(min_val,max_val,n);
  if size(z,1) > 1
    z = z';
  endif
 
  for i = 1:n
    u = (z(i) - x)./h;
    K = (1/sqrt(2*pi))*exp(-.5*u.^2);
    f(i) = sum(y.*K)/sum(K);
  endfor

  if (nargout == 2)
    yy = f;
    xx = z;
  else
    plot(z,f);
  endif

endfunction