nlreg.m

matlab 非线性回归算法 是对现有数据系列进行参数拟合的程序

%
% 2D Nonlinear Regression Fit of y = f(x)
%
%      [ p, resid, h ] = nlreg( y, x, po, pr, 'func' )
%
% where  p      = fitted parameters
%        resid  = residuals
%        h      = handles of the five plots:
%                 h(1) = data and regression line
%                 h(2) = distribution of residuals
%                 h(3) = residuals against x
%                 h(4) = residuals against y
%                 h(5) = parameter confidence region
%        y      = dependent data vector
%        x      = independent data vector
%        po     = initial parameter estimates
%                 ( one set per row )
%        pr     = range for confidence region plot
%                 [ min_p(1) max_p(1) min_p(2) max_p(2) ]
%        'func' = string with model function
%                 ( see below )
%
% the size of y and x should be the same
%
% The following is an example of the model function which
% must put the calculated dependent variable into a global
% array called "yhat" and return the sum-of-squares 
%
%        function sos = monodfn( p )
%        global mudat
%        global sdat
%        global yhat
%        yhat = p(1) .* sdat ./ ( p(2) + sdat );
%        resid = mudat - yhat;
%        sos = resid * resid';
%        return;
%
% Bob Newell, February 1996
%
%------------------------------------------------------
%
function [ p, resid, h ] = nlreg( y, x, po, pr, func )
%
global yhat
%
fontsz = 16;
%
y = y(:);
x = x(:);
disp( 'NONLINEAR REGRESSION' );
h = zeros( 5, 1 );
[ no, np ] = size( po );
%
% optimisation performed here - lots of work
% because we may be doing it from different
% starting points just to be sure we are
% really converging on the minimum
%
options = [ 0 0.1 0.1 ];
p = [];
evalstr = [ 'pnl = fmins( ''' func ''', po(i,:), options );' ];
for i = 1:no,
  eval( evalstr );
  p = [ p; pnl ];
end;
%--------------------------------------------------------
% show the results
disp( 'Optimisation Results:' )
disp( p );
% use the average of the results from here
if no > 1,
  pnl = mean( p );
else,
  pnl = p;
end;
disp( 'Fitted Parameters:' )
for i = 1:np,
  disp( [ '     p(' num2str( i ) ') = ' num2str( pnl(i) ) ] )
end;
% calculate the residuals at the average parameters
% ( results left in global variable )
sos = feval( func, pnl );
% plot the results
h(1) = figure;
plot( x, yhat, x, y, '+' )
set( gca, 'fontsize', fontsz );
xlabel( 'X Data' )
ylabel( 'Y Data' )
%---------------------------------------------------------
% residuals plots
%
yhat = yhat(:);
resid = y - yhat;
%
% plot the distribution first
h(2) = resdist( resid );
set( gca, 'fontsize', fontsz );
xlabel( 'Residuals' )
ylabel( 'Frequency' )
% plot against x data
h(3) = figure;
plot( x, resid, '*' )
set( gca, 'fontsize', fontsz );
ax = axis;
line( [ ax(1:2) ], [ 0 0 ] )
xlabel( 'X Data' )
ylabel( 'Nonlinear Residuals' )
% plot against y data
h(4) = figure;
plot( y, resid, '*' )
set( gca, 'fontsize', fontsz );
xlabel( 'Y Data' )
ylabel( 'Nonlinear Residuals' )
ax = axis;
line( [ ax(1:2) ], [ 0 0 ] )
%-------------------------------------------------
% percentage of variation explained
%
% total variation
sostot = ( y - mean(y) )' * ( y - mean(y) );
pcfit = 100 * ( sostot - sos ) / sostot;
disp( [ 'Percent Variance Explained by Fit = ' ...
             num2str(floor(pcfit)) ] ) ;
%-------------------------------------------------
% parameter confidence region
%
if length( pr ) == 0, return; end;
%
% read from tables, for 95% and dof=(2,96), F=7.53
%
% sos on the 95% contour
sos95 = sos * ( 1 + 2 * 7.53 / 96 );
% plot these contours
crv = [ 1:4 ] * sos95;
% grid of parameter values
crx = [ 1:21 ] * ( pr(2) - pr(1) ) / 20 + pr(1);
cry = [ 1:21 ] * ( pr(4) - pr(3) ) / 20 + pr(3);
% calculate the sos values
s = zeros( 21, 21 );
evalstr = [ 's(i,j) = ' func '( [ crx(i) cry(j) ] );' ];
for i = 1:21,
  for j = 1:21,
    eval( evalstr );
  end;
end;
% plot the contours
h(5) = figure;
contour( crx, cry, s, crv )
set( gca, 'fontsize', fontsz );
grid
xlabel( 'Parameter 1' )
ylabel( 'Parameter 2' )
% add a blob for the parameter estimates
hold on
ax = axis;
blobx = 0.020 * ( ax(2) - ax(1) );
bloby = 0.025 * ( ax(4) - ax(3) );
xblob = blobx * [ 1 -1 -1 1 ]' + pnl(1);
yblob = bloby * [ -1 -1 1 1 ]' + pnl(2);
fill( xblob, yblob, 'k' )
hold off
%
return;
%------------------------------------------------
% the end