<< linregr Regression tools (regtools) nlinlsq(Conflicted) >>

Regression tools (regtools) >> Regression tools (regtools) > nlinlsq

nlinlsq

Non-linear [weighted] least square solver with statistical analysis of the solution.

Calling Sequence

[p,[stat]]=nlinlsq(fun,dfun,y,[wt],p0,[plo],[pup],[info],[algo],[df],[alfa],[pscaling])
p=nlinlsq(fun,[],y,[],p0) // shortest call - with numerical derivatives
nlinlsq(stat)  // print report from a previously solved regression problem.

Parameters

fun:

list - list(fundef,x,[optional parameters])

fundef:

name of non-linear model with calling sequence fundef(p,x,[optional parameters]), where p is a vector of regression parameters.

x:

matrix - independent variables (as column vectors)

dfun:

optional list - list(dfundef,x,[optional parameters] defining analytical derivatives. If dfun=[] then numerical derivatives are calculated.

dfundef:

function calculating the analytical derivative of fundef wrt the regression parameters p.

y:

column vector - dependent variables

wt:

optional column vector with weights for each dependent variable y (typically wt=(1)./y).

p0:

column vector - inital estimates for regression parameters

plo:

optional column vector - lower bounds for regression parameters

pup:

optional column vector - upper bounds for regression parameters

info:

optional list - output options

info(1):

scalar - =0; (default) no output; =1 summary report after solution; >1 output after every info(1)'th iteration.

info(2):

optional text matrix or space separated string with names for the regression parameters. If not present p(1), p(2),... are used as names for the regression parameters in p.

algo:

optional string - solution method ('qn' -(default) quasi-newton, 'gc' - conjugate gradient, 'nd' - non-differentiable model.)

df:

optional scalar - number of degrees of freedom (default df=length(y)-length(p0))

alfa:

optional significance level for parameter confidence interval estimates (default alfa=0.05).

pscaling:

optional scaling factors to improve convergence (default=ones(p0)). Try values close to p.

p:

vector - least square solution of regression parameters p

stat:

optional structure - statistical data from parmeter estimation problem.

stat.ss :

weighted residual sum of squares ( ss = (wt.*res)'*(wt.*res) )

stat.df :

degrees of freedom ( df = length(y) - length(p) )

stat.res :

vector with residuals ( res = y-funlist(p,x,...) )

stat.p :

vector with solution of the WLSQ problem ( min sum(res.^2) subject to plo le p le pup )

stat.pint :

confidence interval (alfa=5%) for reqression parameters ( pint = devp*cdft('T',df,1-alfa/2,alfa/2) )

stat.covp :

parameter covariance matrix ( covp = inv(df/ss*(J'*J)) where J is the Jacobi matrix of res wrt p at the solution)

stat.corp :

parameter correlation matrix ( corp = covp./sqrt(devp2*devp2'); where devp2=devp.^2; )

stat.devp :

standard error ( devp = sqrt(diag(covp)) ).

Description

nlinlsq solves non linear weigthed least square problems using the Scilab function optim as the non-linear optimization routine.

Minimize SS = SUM(i, (wt*(y(i) - f(p,x(i))))^2 ) with respect to plo le p le pup.

The gradient information for the optim algorithm is estimated numerically using the numdiff function, unless dfun is defined - in which case df/dp = dfundef(p,x,...) is used. dfundef(p,x,...) should return [df/dp(1) df/dp(2) ... df/dp(np)] as a column matrix.

Model parameters are automatically scaled using the values in p0 as scaling factors. If max(p/p0) gt 1000 or min(p/p0)le 1e-3 it may be wise to re run nlinlsq with the new solution as initial quess in order to provide better scaling of the model parameters.

Statistical analysis inspired by the nls function in R (www.r.org) is available through the output variable stat.

Examples

deff('yhat=nlinmod(p,x)',['yhat=p(1)+p(2)*exp(p(3)*x);']); // define regression model to be fitted
deff('dydp=dnlinmod(p,x)',['dydp=[ones(x),exp(p(3)*x),p(2)*x.*exp(p(3)*x)];']); // define d(yhat)/dp
x=-[1:100]'/10; phat=[2;10;0.5];               // generate some data points
y=nlinmod(phat,x)+grand(100,1,'nor',0,1)/2;    // with with random noice added.
p0=ones(3,1); // initial estimate for the regression parameters.
// Solve nonlinear regression problem with output every 4'th iteration and nameing of model parmameters.
[p,stat]=nlinlsq(list(nlinmod,x),list(dnlinmod,x),y,[],p0,[],[],list(4,'A B C'));

// Solve weighted nonlinear regression problem with default names for the regression parameters
// and numerical derivatives.
[pwt,stat]=nlinlsq(list(nlinmod,x),list(dnlinmod,x),y,(1)./y,p0,[],[],10);

// Show the difference between the two solutions...
scf(); plot(x,y,'o'); xtitle('Demo of nlinlsq()','x','y=A+B*exp(C*x)')
plot(x,nlinmod(p,x),'b-'); plot(x,nlinmod(pwt,x),'r-');
xgrid(); legend('Data','unweighted','weighted',2);

// Solve weighted nonlinear regression problem without analytical derivaties.
[pwt,stat]=nlinlsq(list(nlinmod,x),[],y,(1)./y,p0,[],[],10);

clc;
// Display the regression report from the previous solution.
nlinlsq(stat)

See also

Authors


Report an issue
<< linregr Regression tools (regtools) nlinlsq(Conflicted) >>