NATL INST. OF STAND 4 TECH

A111DS TbflTSfi

4

V

\

NBS

Publications

*"*e*v of * NBS TECHNICAL NOTE 1068-2

U.S. DEPARTMENT OF COMMERCE / National Bureau of Standards

109

U5753

1063-2

1933
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards' was established by an act ot Congress on March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and technology

and facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau's technical work is per- formed by the National Measurement Laboratory, the National Engineering Laboratory, and

the Institute for Computer Sciences and Technology.

THE NATIONAL MEASUREMENT LABORATORY provides the national system of

physical and chemical and materials measurement; coordinates the system with measurement

systems of other nations and furnishes essential services leading to accurate and uniform

physical and chemical measurement throughout the Nation's scientific community, industry,

and commerce; conducts materials research leading to improved methods of measurement,

standards, and data on the properties of materials needed by industry, commerce, educational

institutions, and Government; provides advisory and research services to other Government

agencies; develops, produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

Absolute Physical Quantities 2 — Radiation Research — Chemical Physics —

Analytical Chemistry — Materials Science

THE NATIONAL ENGINEERING LABORATORY provides technology and technical ser- vices to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement capabilities;

provides engineering measurement traceability services; develops test methods and proposes

engineering standards and code changes; develops and proposes new engineering practices;

and develops and improves mechanisms to transfer results of its research to the ultimate user. The Laboratory consists of the following centers:

Applied Mathematics — Electronics and Electrical Engineering2 — Manufacturing

Engineering — Building Technology — Fire Research — Chemical Engineering2

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts

research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and

economy in Government operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing the Federal Information Processing Standards Program, developing Federal ADP standards

guidelines, and managing Federal participation in ADP voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and

provides the technical foundation for computer-related policies of the Federal Government.

The Institute consists of the following centers:

Programming Science and Technology — Computer Systems Engineering.

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address Washington, DC 20234. 'Some divisions within the center are located at Boulder, CO 80303.
flATIO:. i. . : AU OF STi il .. . DS

LIBi . I

Nonlinear Least Squares Regression ^

Using STARPAC ' 6

^

x^%i

The Standards Times Series

and Regression Package

Janet R. Donaldson

Peter V. Tryon

Statistical Engineering Division

Center for Applied Mathematics

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

c

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Issued October 1 983
National Bureau of Standards Technical Note 1068-2

Natl. Bur. Stand. (U.S.), Tech. Note 1068-2, 64 pages (Oct. 1983)

CODEN: NBTNAE

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1983

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC 20402

Price $4.50

(Add 25 percent for other than U.S. mailing)
In memory of

Peter V. Tryon

1941 - 1982
Disclaimer

No warranties 3 express or implied, are made

by the distributors or developers that

STARPAC or its constituent parts are free of

error. "Vh.ey should not be relied upon as the sole basis for solving a problem whose

incorrect solution could result in injury to

person or property . If the programs are

employed in such a manner, it is at the

user's ovm risk and the distributors and

developers disclaim all liability for such

misuse.

Computers have been identified in this paper

in order to adequately specify the sample

programs and test results. Such

identification does not imply recommendation

or endorsement by the National Bureau of

Standards, nor does it imply that the

equipment identified is necessarily the best

available for the purpose.

iv
Preface

STARPAC, the Standards Time Series and Regression Package, is a library

of Fortran subroutines for statistical data analysis developed by the

Statistical Engineering Division (SED) of the National Bureau of Standards

(NBS), Boulder, Colorado. Earlier versions of this library were distributed

by the SED under the name STATLIB [Tryon and Donaldson, 1978] . STARPAC

incorporates many changes to STATLIB, including additional statistical

techniques, improved algorithms and enhanced portability.

We are distributing STARPAC source code in segments, and the

documentation as a series of NBS Technical Notes. This will allow us to make

the new capabilities of STARPAC available more quickly and to facilitate the

release of future changes.

The first of the series of Technical Note documenting STARPAC is an Introduction to STARPAC [Donaldson and Tryon, 1983] which gives an an overview

of the STARPAC library, defines conventions used to specify the documentation,

and presents general background material. The Introduction to STARPAC incudes

information which is essential for using the STARPAC library, and users should

be familiar with its contents before attempting to use any STARPAC

subroutine.

This Technical Note documents the STARPAC nonlinear least squares

regression subroutines. These subroutines:

• use the quasi-Newton NL2S0L algorithm [Dennis et al. , 1981a] , designed

specifically for problems which may have highly nonlinear models, large

residuals, or poor scaling;

• accept user-supplied derivatives (the Jacobian matrix) when available,

and approximate the derivatives numerically otherwise;

• provide utility procedures to verify the correctness of the

user-supplied derivatives, or to select optimum step sizes for

numerically approximating the derivative [Schnabel, 1982];

• allow both weighted and unweighted analyses;

• permit subsets of the model parameters to be treated as constants with

their values held fixed at their input values, allowing the user to

examine the results obtained by estimating subsets of the parameters of

a general model without rewriting the model subroutine;

• are easy to use, providing three levels of user-control of the

computations and results, extensive error handling facilities,

comprehensive printed reports, and no size restrictions other than

effective machine size;
• are portable; and

• are easily used with other Fortran subroutine libraries.

Other code segments of STARPAC include subroutines for time series

analysis (in both time and frequency domains), line printer graphics, basic

statistical analysis, linear least squares regression, and nonlinear

optimization. Users may obtain a list of the available code segments and

their associated documentation from their local installer of STARPAC, or directly from the author. Notation, format, and naming conventions will

remain constant throughout the series of Technical Notes documenting STARPAC,

allowing the documentation for each code segment to be used alone or in

conjunction with the documentation for another.

The code for STARPAC was developed at the U.S. Department of Commerce

Boulder Laboratories on a CDC CYBER 170/750 under NOS 1.4, and all examples

presented in this documentation were executed in this environment using the

FTN 4.8+552 compiler with rounding. STARPAC has also been tested on the

following equipment.

Computer Facility

IBM 4 340

VAX 11/780

National Center for Atmospheric Research,

Boulder, Colorado

Sperry 1100/82 National Bureau of Standards,

Gaithersburg, Maryland

Perkin-Elmer 3230 National Bureau of Standards,

Boulder, Colorado

Data General ECLIPSE S/140 Forest Fire Laboratory,

Riverside, California

STARPAC is written in Fortran '66. Adherence to a portable subset of

ANSI Fortran [1966] has been verified by the PFORT Verifier [Rvder, 1974].

Workspace and machine-dependent constants are supplied using subroutines based

on the Bell Laboratories "Framework for a Portable Library" [Fox et al.

,

1978a]. We have also used subroutines from UNPACK [Dongarra et al. , 1979],

from the "Basic Linear Algebra Subprograms for Fortran Usage" [Lawson et al.

,

1979], and from DATAPAC [Filliben, 1977]. The subroutines used to compute the

least squares solution are those referenced in Dennis et al. [1981a and

1981b], and the algorithms used to select optimum step sizes for numerical

derivatives, and to check analytic derivatives were developed by Schnabel

[1982]. The printed reports from the nonlinear least squares subroutines have

been adapted from OMNITAB II [Hogben et al. , 1971].

VI
Computer facilities for the STARPAC project have been provided in part by

the NOAA Computer Services Division, Boulder, Colorado, and we gratefully

acknowledge their support. The STARPAC subroutine library is the result of

the programming efforts of Janet R. Donaldson and John E. Koontz, with

assistance from Ginger A. Caldwell, Steven M. Keefer, and Linda L. Mitchell.

Valuable contributions have also been made by each of the members of the

Statistical Engineering Division in Boulder, and from many within the STARPAC

user community. We are grateful for the many valuable comments that we have

received on early drafts of the STARPAC documentation; we wish especially to

thank Paul T. Boggs, Ginger A. Caldwell, Sally E. Howe, John E. Koontz, James

T. Ringland, Robert B. Schnabel, Ralph M. Slutz, and Dominic F. Vecchia.

Finally, we wish to thank Lorna Buhse for excellent manuscript support.

Janet R. Donaldson

Peter V. Tryon (deceased)

September 1983

vii
Contents

Preface v

1 . Nonlinear Least Squares Regression. 2-1 .1

A. Introduction 2-1.1

B. Subroutine Descriptions 2-1.2

B.l Nonlinear Least Squares Estimation Subroutines 2-1.2

B.2 Derivative Step Size Selection Subroutines 2-1.3

B.3 Derivative Checking Subroutines 2-1.4

C. Subroutine Declaration and CALL Statements 2-1.5

D. Dictionary of Subroutine Arguments and

COMMON Variables 2-1.10

E . Computational Methods 2-1.20

E .1 Nonlinear Least Squares Estimation 2-1.20

E.l.a The Nonlinear Least Squares Algorithm 2-1.20

E.l.b Computed Results and Printed Output From

the Nonlinear Least Squares Regression

Subroutines 2-1.22

E .2 Step Size Selection 2-1.25

E.2.a The Step Size Selection Algorithm ."T. . . . . 2-1.25

E.2.b Computed Results and Printed Output From

the Derivative Step Size Selection

Subroutines 2-1.28

E .3 Derivative Checking 2-1.28

E.3.a The Derivative Checking Algorithm 2-1.28

E.3.b Computed Results and Printed Output From

the Derivative Checking Subroutines 2-1.29

F . Examp les 2-1 . 3

Appendix A 2-A.l

Appendix B 2-B .1

References 2-R.l

viii
Nonlinear Least Squares Regression

Using STARPAC

The Standards Time Series and Regression Package

Janet R. Donaldson and Peter V. Tryon

National Bureau of Standards

Washington, DC 20234

STARPAC, the Standards Time Series and Regression Package, is a library of Fortran subroutines for statistical data analysis

developed by the Statistical Engineering Division (SED) of the

National Bureau of Standards (NBS), Boulder, Colorado. Earlier

versions of this library were distributed by the SED under the name STATLIB [Tryon and Donaldson, 1978]. STARPAC incorporates many

changes to STATLIB, including additional statistical techniques,

improved algorithms and enhanced portability.

STARPAC emphasizes the statistical interpretation of results,

and, for this reason, comprehensive printed reports of auxiliary

statistical information, often in graphical form, are automatically

provided to augment the basic statistical computations performed by

each user-callable STARPAC subroutine. STARPAC thus provides the

best features of many stand-alone statistical software programs

within the flexible environment of a subroutine library.

This Note documents 16 subroutines for nonlinear least squares

regression. Twelve of these compute the least squares estimates,

performing either weighted or unweighted analysis with either

numerically approximated or user-supplied (analytic) derivatives.

The other four are user-callable subroutines for two procedures used

within the estimation code: the first selects optimum step sizes

for approximating the partial derviatives of the model; and the

second checks the validity of a user-supplied derivative

subroutine.

Key words: derivative checking; NL2S0L; nonlinear least squares;

nonlinear regression; quasi-Newton methods; STARPAC; statistical

computing; statistical subroutine library; statistics; derivative

step size selection; weighted nonlinear least squares.

ix

CHAPTER 1

NONLINEAR LEAST SQUARES REGRESSION

A. Introduction

STARPAC contains 16 user-callable subroutines for nonlinear least squares

regression. Twelve of these are estimation subroutines that compute the least

squares solution as described below, performing either weighted or unweighted

regression with either numerically approximated or user-supplied (analytic)

derivatives. The estimation subroutines allow three levels of control of the

computations and printed output, and allow the user to specify a subset of the

parameters to be treated as constants, with their values held fixed at their

input values. This last feature allows the user to examine the results

obtained estimating various subsets of the parameters of a general model

without rewriting the model subroutine for each subset. The other four

subroutines described in this chapter are utility procedures which choose

optimum step sizes for numerically approximating the derivative, and which

verify the correctness of user-supplied (analytic) derivatives.

Each of the subroutines described in this chapter assumes that the

observations of the dependent variable, y^ , are modeled by

y±

= f i( x i>£) + £± for 1-1, .-., N,

where

N is the number of observations;

f^ is the function (nonlinear in its parameters) that models the i tn

observation;

x^ is the vector of the M independent variables at the i tn observation;

3 is the vector of the NPAR model parameters; and

e^ is the unobservable random error in the i tn observation, which is

estimated by the i tn residual.

The least squares estimates of the parameters, 6, are obtained using an iterative procedure that requires the matrix of partial derivatives of the

model with respect to each parameter,

D(i,k) = 3f1 (x 1 ,3)/33(k) for 1=1, ..., N

and k = 1, . .

. , NPAR.

The derivative matrix may be supplied analytically or approximated

numerically

.

2-1.1
The least squares solution, 3, is that which minimizes (with respect to

3) the residual sum of squares function,

N N

RSS(3) = I [{yi ~ f

± (x it &)) 2 -wt± ] =

I [e ± 2-wt± l i=l 1=1

where carat (*) denotes the estimated quantity, and

wt^ is the weight assigned to the i

fc " observation (wt^ = 1.0 in the

"unweighted" case). Appendix A discusses several common applications

for weighted least squares.

The user must supply both initial values for the parameters, and a subroutine [see §D, argument NLSMDL] to compute f^x^,^), 1 = 1, . . . , N,

i.e., the predicted values of the dependent variable given the independent

variables and the parameter values from iteration I, % = 1, 2, ... . Initial

parameter values should be chosen with care, since good values can significantly reduce computer time.

STARPAC provides a variety of subroutines to accommodate many levels of

user sophistication and problem difficulty. Users are directed to §B for a brief description of the features and subroutines available. The actual

declaration and CALL statements are given in §C, and the subroutine arguments

are defined in §D. Generally, the computational details provided in §E are not needed unless difficulties arise. Sample programs and their output are

shown in %¥

.

B . Subroutine Descriptions

B .1 Nonlinear Least Squares Estimation Subroutines

The simplest of the 12 nonlinear least squares estimation subroutines,

NLS, requires neither user-supplied weights nor analytic derivatives. The

estimated results and a variety of statistics are automatically summarized in

a five-part printed report, and the estimated parameters and residuals are returned to the user via the subroutine argument list (level one control,

described below). Most nonlinear least squares problems can be solved using

NLS.

The other 11 estimation subroutines add the weighting, derivative, and

level two and three control features both singly and in combination, providing

greater flexibility to the user at the price of less simplicity. These

features are indicated by the suffix letter(s) on the subroutine name (e.g.,

NLSS and NLSWDC).

2-1.2
• Suffix W indicates user-supplied weights.

• Suffix D indicates user-supplied (analytic) derivatives.

• Suffix C indicates level two control of the computations.

m Suffix S indicates level three control of the computations.

The three levels of computation and printed output control are described

as follows.

o In level one, a five-part printed report is automatically provided

[see §E.l.b], and the estimated model parameters and residuals are returned to the user via the argument list.

• Level two also returns the estimated parameters and residuals, and,

in addition, allows the user to supply arguments to indicate

- a subset of the model parameters to be treated as constants,

with their values held fixed at their input values;

- either the step sizes used to compute the numerical

approximations to the derivative [see §E.2], or, when

user-supplied analytic derivatives are used, whether they will

be checked [see §E.3];

- the maximum number of iterations allowed;

- the convergence criteria;

- the scale (i.e., the typical size) of each parameter;

- the maximum change allowed in the parameters at the first

iteration;

- how the variance-covariance matrix is to be approximated; and

- the amount of printed output desired.

• Level three has all the features of level two, and, in addition

returns the following estimated values via the argument list:

- the number of nonzero weighted observations (only when a weighted analysis is performed);

- the number of parameters actually estimated;

- the residual standard deviation;

- the predicted values;

- the standard deviations of the predicted values;

- the standardized residuals; and

- the variance-covariance matrix of the estimated parameters.

B .2 Derivative Step Size Selection Subroutines

When the partial derivatives used in the nonlinear least squares solution

are not available analytically, STARPAC subroutines approximate them

numerically. In this case, the subroutines can select optimum step sizes for

approximating the derivatives [see §E.2]. The user also has the option of

computing these step sizes independently of the estimation process by calling

either of the two step size selection subroutines directly. For example, when

2-1.3
planning to use the parameter fixing capability [see §D, argument IFIXED] to

examine several subsets of the parameters of a general model, computing the

step sizes first and passing them to the estimation subroutine is more

efficient than recomputing them each time the estimation subroutine is

called.

The simplest of the two user-callable step size selection subroutines,

STPLS, summarizes the step size selection information for each parameter in a printed report, and returns the step sizes to the user via the subroutine

argument list.

The second step size selection subroutine, STPLSC, differs from STPLS

only in that it enables the user to supply arguments to indicate

- the number of reliable digits in the model results;

- the number of exemptions allowed by the acceptance criteria;

- the scale (i.e., the typical size) of each parameter; and

- the amount of printed output desired.

B .3 Derivative Checking Subroutines

When the partial derivatives used in the nonlinear least squares solution

are available analytically, the user can code them for use by the estimation

subroutines [see §D, argument NLSDRV] . Because coding errors are a common

problem with user-supplied derivatives, the STARPAC estimation subroutines

automatically check the validity of the user-supplied derivative code by

comparing its results to numerically approximated values for the derivative.

When the results are questionable, the checking procedure attempts to

determine whether the problem lies with the user's code or with the accuracy

of the numerical approximation [see §E.3]. Although the checking procedure is

automatically available to the estimation subroutines which accept

user-supplied derivatives, the user may want to check the derivative code

independently of the estimation process. In these cases, the user can call

either of the two derivative checking subroutines directly, and suppress

checking by the estimation subroutines. [See §D, argument IDRVCK.

]

The simplest of the two derivative checking subroutines, DCKLS,

summarizes the results of the check in a printed report.

The second derivative checking subroutine, DCKLSC, differs from DCKLS

only in that it enables the user to supply arguments to indicate

- the number of reliable digits in the model results;

- the agreement tolerance;

- the scale (i.e. , the typical size) of each parameter;

- the row at which the derivative is to be checked; and

- the amount of printed output desired.

2-1.4
C . Subroutine Declaration and CALL Statements

NOTE: Argument definitions and sample programs are given in §D and §F

,

respectively. The conventions used to present the following declaration and

CALL statments are given in the Introduction to STARPAC [Donaldson and Tryon,

1983, chapter 1, §B, and §D] .

Nonlinear Least Squares Estimation Subroutines

The <basic declaration block> identifies declaration statements that are needed by all of the nonlinear least squares estimation subroutines. The user should substitute the following four statements for each occurrence of <basic

declaration block> given below.

<real> Y(n), XM(n,r?), ?AR(npar) , RES(^?)

DOUBLE PRECISION DSTAK( Idstak)

COMMON /CSTAK/ DSTAK

EXTERNAL NLSMDL

NLS: Unweighted nonlinear least squares regression; numerically

approximated derivatives; five-part printed report; minimal return of

results.

<basic declaration block>

CALL NLS (Y, XM , N, M, IXM , NLSMDL, PAR, NPAR, RES, LDSTAK)

NLSC: Unweighted nonlinear least squares regression; numerically

approximated derivatives; user-supplied control values; minimal return

of results.

<basic declaration block>

INTEGER IFIXED( npar)

<real> ST? (n par) , STOPSS, STOPP, SCALE(npar) , DELTA

CALL NLSC (Y, XM , N, M, IXM, NLSMDL, PAR, NPAR, RES, LDSTAK,

1 IFIXED, STP, MIT, STOPSS, STOPP, SCALE, DELTA, IVAPRX, NPRT)

2-1.5
NLSS: Unweighted nonlinear least squares regression; numerically

approximated derivatives; user-supplied control values; maximal return

of results.

<basic declaration block>

INTEGER IFIXEV( npar)

<real> ST? (n par) , STOPSS, STOPP, SCKLE(npar) , DELTA

<real> RSD, PV( n) , SDPV( n) , SDRES(w), VCV(npare3 npare)

CALL NLSS (Y, XM , N, M, IXM , NLSMDL, PAR, NPAR, RES, LDSTAK,

1 IFIXED, STP, MIT, STOPSS, STOPP, SCALE, DELTA, IVAPRX, NPRT,

2 NPARE, RSD, PV, SDPV, SDRES, VCV, IVCV)

NLSW: Weighted nonlinear least squares regression; numerically approximated

derivatives; five-part printed report; minimal return of results.

<basic declaration block>

<real> WT(rz)

CALL NLSW (Y, WT, XM , N, M, IXM, NLSMDL, PAR, NPAR, RES, LDSTAK)

NLSWC: Weighted nonlinear least squares regression; numerically approximated

derivatives; user-supplied, control values; minimal return of results.

<basic declaration block>

INTEGER IFIXED(rzpar)

<real> WT(rz)

<real> STP (n par) , STOPSS, STOPP, SCALE( npar) , DELTA

CALL NLSWC (Y, WT, XM , N, M, IXM, NLSMDL, PAR, NPAR, RES, LDSTAK,

1 IFIXED, STP, MIT, STOPSS, STOPP, SCALE, DELTA, IVAPRX, NPRT)

2-1.6
NLSWS: Weighted nonlinear least squares regression; numerically approximated

derivatives; user-supplied control values; maximal return of results,

<basic declaration block>

INTEGER IFIXED( n par)

<real> WT(n)

<real> STP(npar), STOPSS, STOPP, SCALE(npar) , DELTA

<real> RSD, PV( n) , SDPV( n) , SDRES(n), VCV ( n par

e

s n pare)

CALL NLSWS (Y, WT, XM , N, M, IXM , NLSMDL, PAR, NPAR, RES, LDSTAK,

1 IFIXED, STP, MIT, STOPSS, STOPP, SCALE, DELTA, IVAPRX, NPRT,

2 NNZW, NPARE, RSD, PV, SDPV, SDRES, VCV, I VCV)

NLSD: Unweighted nonlinear least squares regression; user-supplied

derivatives; five-part printed report; minimal return of results,

<basic declaration block>

EXTERNAL NLSDRV

CALL NLSD (Y, XM, N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK)

NLSDC: Unweighted nonlinear least squares regression; user-supplied

derivatives; user-supplied control values; minimal return of results,

<basic declaration block>

EXTERNAL NLSDRV

INTEGER IFIXED( npar)

<real> STOPSS, STOPP, SCALE(npar), DELTA

CALL NLSDC (Y, XM , N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK, IFIXED, IDRVCK, MIT, STOPSS, STOPP, SCALE, DELTA,

2 IVAPRX, NPRT)

2-1.7
NLSDS: Unweighted nonlinear least squares regression; user-supplied

derivatives; user-supplied control values; maximal return of results,

<basic declaration block>

EXTERNAL NLSDRV

INTEGER IFIXED(npa^)

<real> STOPSS, STOPP, SCALE(npar), DELTA

<real> RSD, PV(w), SDPV( n) , SDRES(n), VCV

n par

3 n pare)

CALL NLSDS (Y, XM , N, M, IXM , NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK, IFIXED, TDRVCK, MIT, STOPSS, STOPP, SCALE, DELTA,

2 IVAPRX, NPRT, NPARE, RSD, PV, SDPV, SDRES, VCV, I VCV)

NLSUD: Weighted nonlinear least squares regression; user-supplied

derivatives; five-part printed report; minimal return of results.

<basic declaration block>

EXTERNAL NLSDRV

<real> WT(rz)

CALL NLSWD (Y, WT, XM , N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK)

NLSWDC: Weighted nonlinear least squares regression; user-supplied

derivatives; user-supplied control values; minimal return of results,

<basic declaration block>

EXTERNAL NLSDRV

INTEGER IFIXED( npar)

<real> WTO)

<real> STOPSS, STOPP, SCALE(npar), DELTA

CALL NLSWDC (Y, WT, XM , N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK, IFIXED, IDRVCK, MIT, STOPSS, STOPP, SCALE, DELTA,

2 IVAPRX, NPRT)

2-1.8
NLSWDS: Weighted nonlinear least squares regression; user-supplied

derivatives; user-supplied control values; maximal return of results,

<basic declaration block>

EXTERNAL NLSDRV

INTEGER IFIXED(npar)

<real> WT(w)

<real> STOPSS, STOPP, SCALE( npar) , DELTA

<real> RSD, PV( n) , SDPV( n) , SDRES(n), VCV (n pare} n pare)

CALL NLSWDS (Y, WT, XM , N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, RES,

1 LDSTAK, IFIXED, IDRVCK, MIT, STOPSS, STOPP, SCALE, DELTA,

2 IVAPRX, NPRT, NNZW, NPARE, RSD, PV, SDPV, SDRES, VCV, IVCV)

Step Size Selection Subroutines

STPLS: Optimum step size selection for numerically approximating

derivatives; printed report; maximal return of results.

<real> XM(n,w), FkR(npar) , ST?(npar)

DOUBLE PRECISION DSTAK( Idstak)

COMMON /CSTAK/ DSTAK

EXTERNAL NLSMDL

CALL STPLS (XM, N, M, IXM, NLSMDL, PAR, NPAR, LDSTAK, STP)

S^PLSC: Optimum step size selection for numerically approximating

derivatives; user-supplied control values; maximal return of results,

<real> m(n,m) , PAR( npar) y STP(npar)

<real> EXMPT, SCALE(npar)

DOUBLE PRECISION DSTAK( Idstak)

COMMON /CSTAK/ DSTAK

EXTERNAL NLSMDL

CALL STPLSC (XM, N, M, IXM, NLSMDL, PAR, NPAR, LDSTAK, STP,

1 NETA, EXMPT, SCALE, NPRT)

2-1.9
Derivative Checking Subroutines

DCKLS: Derivative checking; -printed report; maximal return of results,

<real> XM(n,m), PAR( npar) t ST?(npar)

DOUBLE PRECISION DSTAK( Idstak)

COMMON /CSTAK/ DSTAK

EXTERNAL NLSMDL, NLSDRV

CALL DCKLS (XM, N, M, IXM , NLSMDL, NLSDRV, PAR, NPAR, LDSTAK)

DCKLSC: Derivative checking; user-supplied control values; maximal return of

results.

<real> XM(n,m), PAR( npar), STP(npar)

<real> SCALE(npar)

DOUBLE PRECISION DSTAK( Idstak)

COMMON /CSTAK/ DSTAK

EXTERNAL NLSMDL, NLSDRV

CALL DCKLSC (XM , N, M, IXM, NLSMDL, NLSDRV, PAR, NPAR, LDSTAK,

1 NETA, NTAU, SCALE, NROW , NPRT)

D. Dictionary of Subroutine Arguments and COMMON Variables

NOTE: —> indicates that the argument is input to the subroutine, and that

the input value is preserved;

< — indicates that the argument is returned by the subroutine; <-> indicates that the argument is input to the subroutine and that

the input value will be overwritten by the subroutine;

indicates that the argument is input to some subroutines and is

returned by others; *** indicates that the argument is a subroutine name;

••• indicates that the variable is passed via COMMON.

<— The matrix of exact dimension N by NPAR that contains the partial

derivatives of the model with respect to each parameter,

PAR(k) , k = 1, ..., NPAR. This argument is used within derivative

subroutine NLSDRV [see argument NLSDRV below].

2-1.10
DELTA —> The maximum scaled change allowed In the parameters at the first

iteration, i.e., 6q [see §E.l.a]. The default value is 100.0.

When DELTA < 0.0, or when DELTA is not an argument of the

subroutine CALL statement, the default value is used. A smaller

value of DELTA may be appropriate if, at the first iteration, the

computation of the predicted values from the user's model

subroutine produces an arithmetric overflow, or the parameters

leave the region of interest in parameter space. A reasonable

alternative to the default value of DELTA is an upper bound to the

scaled change that the estimated parameters should be allowed to

make on the first iteration,

DELTA = min{ |AmaxB(k) |/SCALE(k), for k = 1, ..., NPAR}

where Amax3(k) is the maximum change allowed for the k

tn parameter

at the first iteration.

DSTAK The DOUBLE PRECISION vector in COMMON /CSTAK/ of dimension at

least LDSTAK. DSTAK provides workspace for the computations. The

first LDSTAK locations of DSTAK will be overwritten during

subroutine execution.

EXMPT —> The proportion used to compute the number of observations,

a = EXMPT»N, for which the forward difference quotient derivative

with respect to a given parameter is exempted from meeting the

acceptance criteria for step size selection [see §E.2.a]. The

default value for EXMPT is 0.1 (10 percent). When the

user-supplied value is outside the range [0.0, 1.0], or when EXMPT

is not an argument of the subroutine CALL statement, the default

value is used.

IDRVCK —> The indicator variable used to designate whether or not the

user-supplied derivative subroutine is to be checked. When

IDRVCK * the derivative is checked, and when IDRVCK =0 it is

not. The default value is IDRVCK * 0. When IDRVCK is not an argument of the subroutine CALL statement, the default value is

used.

IERR An error flag returned in COMMON /ERRCHK/. [See Donaldson and

Tryon, 1983, chapter 1, §D.5.] Note that using (or not using) the

error flag will not affect any error messages that are automatically printed, even when the user has suppressed the

normal printed output.

For the estimation subroutines:

IERR = indicates that no errors were detected, and that the

iterations converged satisfactorily.

IERR = 1 indicates that improper input was detected.

continued —

2-1.11
IERR = 2 indicates that the computation of the residual sum of

squares using the initial parameter values produced an arithmetic overflow. The user should reduce the size

of DELTA or should supply new starting values.

IERR = 3 indicates that the model is computationally singular,

which means the model has too many parameters near the

solution. The user should examine the model and data

to determine and remove the cause of the singularity.

IERR = 4 indicates that at least one of the standardized

residuals could not be computed because its standard

deviation was zero. The validity of the covariance

matrix is questionable.

IERR = 5 indicates false convergence. [See §E.l.a.]

IERR = 6 indicates that convergence was not reached in the

allowed number of iterations or model subroutine

calls. [See argument MIT.]

IERR = 7 indicates that the variance-covariance matrix could

not be computed.

For the step size selection subroutines:

IERR = indicates that no errors were detected, and that all

the step sizes satisfied the selection criteria.

IERR = 1 indicates that improper input was detected.

IERR = 2 indicates that one or more of the step sizes did not

satisfy the selection criteria.

For the derivative checking subroutines:

IERR = indicates that no errors were detected, and that the

user-supplied derivative code appears to be correct.

IERR = 1 indicates that improper input was detected.

IERR = 2 indicates that the user-supplied derivative code and

numerical derivatives do not agree for at least one

parameter, but that in each case of disagreement the

accuracy of the numerical derivatives is questionable.

Further testing is suggested.

IERR = 3 indicates that the user-supplied derivative code and

numerical derivatives do not agree for at least one

parameter, and in at least one instance of

disagreement there is no reason to doubt the numerical

derivatives.

IFIXED —> The vector of dimension at least NPAR that contains values used to

indicate whether the corresponding parameter in PAR is to be

treated as a fixed constant or is to be estimated. If

IFIXED(I) > 0, PAR(I) will be held fixed at its input value; if

IFIXED(I) = 0, PAR(I) will be estimated using the least squares

continued —

2-1.12
IVCV

procedure described in §A. The default values are IFIXED(I) = 0,

I = 1, •••, NPAR, i.e., all parameters are estimated. When

IFIXED(l) < -1, or when IFIXED is not an argument of the

subroutine CALL statement, the default value will be used.

•> The exact value of the first dimension of the matrix VCV as specified in the calling program.

IVAPRX —> The indicator variable used to specify how the variance-covariance

matrix, VCV, is to be approximated. Three approximations are available:

(1) VCV = RSD 2 «(DT .W.D)

(2) VCV = RSD 2 .H~1

(3) VCV = RSD 2 .H_1 .(DT .W.D).H-1

where

H = D

T .W.D + { y ei'wt. .(8 2

e i /3B(j)96(k)) for j = 1, ..., NPAR

i=l and k = 1, . .

.

, NPAR}

W is an N by N diagonal matrix of user-supplied weights,

W = diag{wtis i = 1, ..., N},

when a weighted analysis is performed, and is the identity matrix

otherwise.

Approximation (1) is based on the assumption that H ~ D »W«D

because the residuals are sufficiently small at the solution;

approximation (2) is based on the assumption that the necessary

conditions for asymptotic maximum likelihood theory have been met;

and approximation (3) is based on the assumption that the

necessary conditions for asymptotic maximum likelihood theory may

be violated. All approximations to the variance-covariance matrix

are subject to sampling variation, however, because they are computed using the estimated parameter values. The

variance-covariance matrix computed for any particular nonlinear

least squares solution should thus be regarded as only a rough

estimate, correct to within an order of magnitude [Bard, 1974].

If IVAPRX = 1 or 4 then approximation (1) is used;

= 2 or 5 then approximation (2) is used; and

= 3 or 6 then approximation (3) is used.

If IVAPRX =1, 2, or 3, then, when user-supplied analytic

derivatives are available [see argument NLSDRV] , they are used to

compute VCV; if IVAPRX =4, 5, or 6, then only the predicted

values from the model subroutine are used to compute VCV. When

analytic derivatives are available, options 1, 2, or 3, will

generally result in a faster, more accurate computation of VCV.

continued —

2-1.13
The default value for IVAPRX is 1. When argument IVAPRX is

outside the range [1, 6], or when IVAPRX is not an argument of the

subroutine CALL statement, then the default value will be used.

IXM —> The exact value of the first dimension of the matrix XM as specified in the calling program.

LDSTAK —> The length of the DOUBLE PRECISION workspace vector DSTAK. LDSTAK

must equal or exceed the appropriate value given below. If the

single precision version of STARPAC is being used, P = 0.5

otherwise P = 1.0 [see Donaldson and Tryon, 1983, chapter 1, $B]

For NLS, NLSC, NLSS, NLSW, NLSWC, NLSWS:

LDSTAK > 27 + max) IS* (N+NPAR) , 30+NPAREJ +

P.max{lS'10.N, 94+N« (3+NPAR)+(3 .NPARE2 +37.NPARE) /2

}

where IS = 1 if default values are used for the derivative

step sizes, and IS = otherwise.

For NLSD, NLSDC , NLSDS, NLSWD, NLSWDC , NLSWDS:

LDSTAK > 45 + NPAR + P« (94+N • (3+NPAR)+(3 .NPARE 2 +35» NPARE) /2 ) For STPLS, STPLSC:

LDSTAK > 27 + (N+NPAR) + P«10-N

For DCKLS, DCKLSC:

LDSTAK > 14 + NPAR + P- (N-NPAR+N+NPAR)

M —> The number of independent variables, i.e., the number of columns

of data in XM.

MIT —> The maximum number of iterations allowed. This argument is also

used to compute the maximum number of model subroutine calls,

(2»MIT). The iterations will stop if either limit is reached,

although, as a rule, the maximum number of iterations will be

reached first. The default value for the maximum number of

iterations is 21. When MIT < 0, or when MIT is not an argument of

the subroutine CALL statement, the default value will be used.

N —> The number of observations.

2-1.14
NETA —> The number of reliable decimal digits in the predicted values (PV)

computed by the user's model subroutine. The default value for

NETA is experimentally determined by the procedure described in

Appendix B. The default value will be used when NETA is not an argument in the subroutine CALL statement, or when the

user-supplied value of NETA is outside the range [1, DIGITS],

where DIGITS is the number of decimal digits carried by the user's

computer for a single precision value when the single precision

version of STARPAC is being used and is the number carried for a double precision value otherwise.

NLSDRV *** The name of the user-supplied subroutine that computes the partial

derivative matrix (Jacobian) . This argument must be listed in an EXTERNAL statement in the program which calls the STARPAC

estimation or derivative checking subroutine. The form of the

derivative subroutine argument list and dimensioning statements

must be exactly as shown below, although if there is only one independent variable (M = 1), XM may be declared to be a vector

with dimension IXM

.

SUBROUTINE NLSDRV (PAR, NPAR, XM , N, M, IXM, D)

<real> PAR(NPAR), XM(IXM,M), D(N,NPAR)

< Computations for D(I,J),

1=1, . . ., N and J = 1 , . .

. , NPAR >

RETURN

END

NLSMDL *** The name of the user-supplied subroutine that computes the

predicted value of the dependent variable given the independent

variables and the current values of the model parameters. This

argument must be listed in an EXTERNAL statement in the program

which calls the STARPAC estimation, step size selection, and/or

derivative checking subroutines. The form of the model subroutine

argument list and dimensioning statements must be exactly as shown

below, although if there is only one independent variable (M = 1)

,

XM may be declared to be a vector with dimension IXM.

SUBROUTINE NLSMDL (PAR, NPAR, XM , N, M, IXM, PV)

<real> PAR(NPAR) , XM(IXM,M), PV(N)

< Computations for PV(I), 1=1, ..., N >

RETURN

END

NNZW <— The number of observations with nonzero weights.

2-1.15
NPAR —> The number of parameters in the model, including both those held

fixed at their starting values, and those which are to be

estimated.

NPARE <— The number of parameters actually estimated, i.e., the number of

zero elements in IFIXED. N.B. This value is returned by the

estimation subroutines.

NPRT —> The argument controlling printed output.

For the estimation subroutines:

NPRT is a five-digit integer, in which the value of the 1^

digit (counting from left to right) is used to control the I tn

section of the output.

If the I tn digit = the output from the I tn section is

suppressed;

= 1 the brief form of the I tn section is

given;

> 2 the full form of the I tn section is

given.

The default value for NPRT is 11112. When NPRT < -1, or when

NPRT is not an argument in the subroutine CALL statement, the

default value will be used. If the convergence criteria are not satisfied, the subroutine gives a suitable warning and

provides a printed report even if NPRT =0. A full discussion

of the printed output is given in §E.l.b, and is summarized as follows.

Section 1 lists the starting estimates and control values.

Brief output and full output are the same for this

section.

Section 2 reports the results of the iterations. Brief output

includes information only about the first and last

iteration while full output includes information

about all of the iterations.

Section 3 provides information for each observation based on

the final solution. Brief output includes

information for the first 40 observations, while full

output provides the information for all of the data.

Section 4 is a set of four residual plots. Brief output and

full output are the same for this section.

Section 5 is the final summary of the estimated parameters.

Brief output does not include printing the

variance-covariance matrix, while full output does.

For the step size selection and derivative checking subroutines:

If NPRT = the output is suppressed.

If NPRT * the output is provided.

If the acceptance criteria were not met, then a printed report

is provided even if NPRT = 0.

2-1.16
NROW —> The row of the independent variable matrix at which the

user-supplied derivative code is to be checked. The default value

is the first row with no independent variables equal to zero; when

all rows have one or more independent variables equal to zero, row one will be used for the default value. When the user-supplied

value is outside the range [1, N] , or when NROW is not an argument

of the subroutine CALL statement, the default value will be used.

NTAU —> The agreement tolerance, i.e., the number of digits of agreement

required between the user-supplied derivatives and the derivatives

numerically approximated by the derivative checking subroutine.

The default value is NETA/4 . When the user-supplied value of NTAU

is outside the range [1, NETA/2], or when NTAU is not an argument

of the subroutine CALL statement, the default value will be" tised.

PAR The vector of dimension at least NPAR that contains the parameter

values. For all estimation subroutines it must contain initial

values for the parameters on input and will contain the final

values on return. For the step size and derivative checking

subroutines it must contain the parameter values at which the

operations are to be performed.

PV <— The vector of dimension at least N that contains the predicted

values of the dependent variable at the solution,

PV(i) = fiCxi.3) = y ±

, for i = 1 , ..

, N.

RES <— The vector of dimension at least N that contains the residuals at

the solution,

RES(i) = y(i) - fiCxj.B) - y(i) - y(i) = e(i) , for i = 1, ..., N.

RSD <— The residual standard deviation at the solution,

RSD = (RSS(3)/(NNZW-NPARE))

1/2

.

SCALE —> The vector of dimension at least NPAR that contains the scale, or typical size, of each parameter. The vector SCALE is used to

normalize the size of each parameter so that

3(j)/SCALE(j) B(k)/SCALE(k) for k = 1, . . ., NPAR

and .1 = 1, • • . , NPAR.

Values of |SCALE(k)| > I B(k) I can be used to increase the step

size in cases where the model function is known to be insensitive

to small changes in the value (3(k) .

For the estimation subroutines:

The default values for SCALE are selected by the NL2S0L

algorithm [Dennis et al. , 1981a and 1981b] and are updated at

each iteration. When SCALE is not an argument in the

2-1.17
subroutine CALL statement, or when the user-supplied value for

SCALE(l) < 0, the default procedure will be used to select

scale values. When SCALE(l) > 0, values of SCALE(k) < 0,

k = 2, ..., NPAR, will be interpreted as an input error. User-supplied scale values may be either a vector of the

typical size of each parameter or a vector of ones if the

typical sizes of the parameters are roughly equal;

user-supplied scale values can sometimes result in reduced

computing time since these values are not updated at each

iteration.

For the derivative checking and step size selection subroutines:

The default values of SCALE are defined as:

SCALE(k) = 1.0 if 3(k) = 0.0

SCALE(k) = l3(k)| otherwise

for k = 1, . .. , NPAR.

When SCALE is not an argument in the subroutine CALL statement, or when the user-supplied value of |SCALE(k)| < l3(k)|, the

default value for SCALE(k) is used. When SCALE(l) < 0, the

default values will be used for SCALE(k), k = 1, ..., NPAR.

When SCALE(l) > 0, values of SCALE(k)

will be interpreted as an input error.

0, k = 2, NPAR,

SDPV

SDRES

<— The vector of dimension at least N that contains an approximation

to the standard deviation of each predicted value at the

solution,

SDPV(i) = the i th diagonal element of [rSD 2 . (D. VCV.Dr

)

1/2

for i = 1, ..., N. This approximation is based on a linearization

of the model in the neighborhood of the solution; the validity of

the approximation depends on the nonlinearity of the model. This

approximation may be extremely inaccurate for a problem with a highly nonlinear model.

<-- The vector of dimension at least N that contains an approximation

to the standardized residuals at the solution,

SDRES(i) = RES(i)/[(RSD 2 /WT(i)) - SDPV(i) 2

] 1/2

for i = 1, . . . , N, which is the i tn residual divided by its

individual estimated standard deviation. This approximation is

based on a linearization of the model in the neighborhood of the

solution; the validity of the approximation depends on the

nonlinearity of the model. This approximation may be extremely

inaccurate for a problem with a highly nonlinear model.

2-1.18
STOPP —> The stopping value for the convergence test based on the maximum

scaled relative change in the parameters at the most recent

iteration. The convergence criterion is satisfied if the current

step is a Newton step [see Dennis et al. , 1981a] and

max{|g8.1 . 1 (k)-3g(k)|/SCALE(k) fork = 1, ..., NPAR} < ST0PP

max{(le £+1 (k)|+|B £ (k)| J/SCALE(k) for k = 1, .... NPAR}

This convergence test is roughly equivalent to the test based on

the maximum relative change in each parameter as measured by

max{ l3£+i(k)-6 £ (k)|/l3 £ (k)| fork = 1, ..., NPAR} . STOPP is not a scale-dependent value; if its value is 10-t+

, then this criteria

will he met when the first four digits of each parameter are the same at two successive iterations, regardless of the size of the

parameter values.

The default value is approximately 10~DIGITS / 2

, where DIGITS is

the number of decimal digits carried by the user's computer for a single precision value when the single precision version of

STARPAC is being used, and is the number carried for a double

precision value otherwise. When the user-supplied value for STOPP

is outside the interval [0.0, 1.0], or when STOPP is not an argument of the subroutine CALL statement, the default value will

be used.

STOPSS —> The stopping value for the convergence test based on the ratio of

the forecasted change in the residual sum of squares, As 2 £+l» to

the current residual sum of squares, s

2 % = RSSCB^) /(NNZW-NPARE)

.

The convergence criterion is satisfied if certain conditions are met [see Dennis et al. , 1981a] and

As 2 £+1 /s 2 £ < STOPSS.

This convergence test is roughly equivalent to the test based on

the relative change in the residual standard deviation between two

iterations % and £+1 as measured by (s^ - s^+i) /$ %• STOPSS is not

a scale-dependent value; if its value is 10~ 5

, this criteria will

be met when the first five digits of the residual sum of squares

are the same at two successive iterations regardless of the size

of the residual sum of squares.

The default value is approximately the maximum of 10 and

10-2 * DIGITS/3 ^ Where DIGITS is the number of decimal digits

carried by the user's computer for a single precision value when

the single precision version of STARPAC is being used, and is the

number carried for a double precision value otherwise. When the

user-supplied value for STOPSS is outside the interval

[10"DIGiTS

, 0.1], or when STOPSS is not an argument of the

subroutine CALL statement, the default value will be used.

2-1.19
STP The vector of dimension at least NPAR that contains the relative

step sizes used to approximate the derivative matrix numerically.

It is input to the estimation subroutines and returned from the

step size selection subroutines. The procedure used to select the

default values is described in §E.2. For the estimation

subroutines, when STP is not an argument of the subroutine CALL

statement or when STP(l) < 0, the default values will be used for

all of the step sizes, and when STP(l) > 0, values of STP(k) < 0,

k = 2, ..., NPAR, will be interpreted as an input error.

VCV <— The matrix of dimension at least NPARE by NPARE that contains the

variance-covariance matrix of the estimated parameters,

approximated as designated by argument IVAPRX. The parameters

which are held fixed [see argument IFIXED] are not included in the

variance-covariance matrix.

The approximation of the variance-covariance matrix is based on a linearization of the model in the neighborhood of the solution;

the validity of the approximation depends on the nonlinearity of

the model. This approximation may be extremely inaccurate for a problem with a highly nonlinear model.

WT —> The vector of dimension at least N that contains the weights.

Negative weights are not allowed and the number of nonzero weights

must equal or exceed the number of parameters being estimated. A

zero weight eliminates the corresponding observation from the

analysis, although the predicted values and standard deviations of

the predicted values are still computed. [See Appendix A.]

XM —> The matrix of dimension at least N by M whose j column contains

the N values of the j

tn independent variable, j = 1, ..., M.

Y —> The vector of dimension at least N that contains the dependent

variable.

E . Computational Methods

E .1 Nonlinear Least Squares Estimation

E.l.a The Nonlinear Least Squares Algorithm

The nonlinear least squares estimation subroutines use the NL2S0L

software package written by Dennis et al. [1981a and 1981b] . The observations

of the dependent variable, which are measured with error, are iteratively fit

to a nonlinear model by minimizing the sums of squares of the errors as described in §A. The iterations continue until the convergence criteria based

on the change in the parameter values or in the residual sum of squares are satisfied [see §D, arguments STOPP and STOPSS] , the maximum number of

iterations (or model subroutine calls) is reached [see §D, argument MIT], or the iterations are terminated due to singularity in the model or false

convergence. All but the first of these stopping conditions may indicate

computational problems and will produce an error report. [See Donaldson and

Tryon, 1983, §D.5.]

2-1.20
Singular convergence means that the model contains too many parameters,

at least near the solution, while false convergence can indicate that either

STOPSS or STOPP is set too small for the accuracy to which the model and its

derivatives are being computed, or that there is an error or discontinuity in

the derivative. Users should examine their models to determine and correct

the underlying cause of singular or false convergence.

Iterative procedures for solving nonlinear least squares problems are discussed in Dennis and Schnabel [1983] , Draper and Smith [1981] , and Kennedy

and Gentle [1980] . The specific procedure used in STARPAC is as follows. At

iteration £+1 the values of the parameter vector $£+i are given by

H+l = H ~ (V-W'D * + S

i + AjJ

_1 -D£T - W ' e

)t T

subject to the restriction that

NPAR

{ I [(0£+l( k > " Mk))/SCALE(k)]2}l/2 < 6£>

k=l

where

NPAR is the number of parameters.

is the vector of NPAR estimated parameter values from the I

i teration.

th

is the N by NPAR matrix of the partial derivatives,

D £ (i,k) = 8f i (x i ,B £ )/36 )l (k) for i = 1, . . ., N

and k = 1, . .

. , NPAR.

W is an N by N diagonal matrix of user-supplied weights,

W = diag{wt^, i = 1, ..., N}, when a weighted analysis is performed,

and is the identity matrix otherwise.

is an approximation to the exact term S^* in the matrix of second

order terms (Hessian) of the Taylor series expansion of the residual

sum of squares function,

N

S **(j> k) =

I

[e £i -wti .(82e £ (i)/3B Ji (j)8B £ (k))],

i=l

e £ SCALE

for j = 1, ..., NPAR, and k = 1, ..., NPAR.

is the vector of the N residuals from the previous iteration.

is the vector of the scale value (i.e., the typical size) of each of

the NPAR parameters. They are either selected by STARPAC or are supplied by the user via subroutine argument SCALE [see §D]

2-1.21
6« is the adaptively chosen diameter of the trust region, i.e., the

region in which the local approximation to the user's model function

is reliable. At each iteration i, 5^ is computed based on information

from the previous iteration. At the first iteration, the initial

value, 6q, is supplied by argument DELTA [see §D] , which can be used

to control the change in the parameters permitted at the first

i teration.

A £ is an NPAR by NPAR diagonal matrix, A

£ = diagj X^/SCALECk) , k = 1, ..., NPAR} , where Ag, is chosen to approximate the smallest non-negative

number such that the restriction given above on the size of the change

in the parameters is satisfied.

The second order term S^ , which is expensive and difficult to compute

accurately, is important only if it is large compared to the term D^^.y.n.,

that is, when the residuals are large or the model is highly nonlinear. When

S^ is large compared to D£T »W»D£, algorithms which ignore it, such as Levenberg-Marquardt or Gauss-Newton, may converge slowly. The NL2S0L

algorithm used by STARPAC, however, adaptively decides when inclusion of this

terra is necessary for reliable results and uses an inexpensive approximation

to Sj£ in those cases.

The matrix, D, of partial derivatives of the model with respect to each

parameter is either computed analytically using a user-supplied subroutine

[see §D, argument NLSDRV], or is numerically approximated using forward

difference quotients as described in §E.2. When the derivatives are approximated numerically, the least squares solution, especially the

variance-covariance matrix, can be sensitive to the step sizes used for the

approximation. The user may want to use STARPAC subroutines STPLS or STPLSC

to recompute the step sizes at the solution provided by the estimation

subroutines to assure that the step sizes which were used are still

acceptable. If there is a significant change in the step size, the least

squares solution should be recomputed with the new step sizes from the current

point. In addition, if the estimation subroutine has convergence problems,

the user may want to recompute the step sizes with the most recent parameter

values to see if a change in the curvature of the model, which will be

reflected as a change in the optimum step sizes, is causing the problem.

Dennis et al. [1981a] provides a detailed description of the algorithm

used in STARPAC. STARPAC also includes the subroutines NL2S0L, NL2SN0, and

NL2ITR, which they reference, and which can be used as documented by them.

[See Dennis et al. , 1981b.]

E.l.b Computed Results and Printed Output

From the Nonlinear Least Squares Regression Subroutines

The argument controlling the printed output, NPRT, is discussed in §D.

The output from the nonlinear least squares estimation subroutines consists of

five sections, several of which include tables summarizing the results. In

the following descriptions, the actual table headings are given by the

underlined, uppercase phrases. Results which correspond to input or returned

subroutine CALL statement arguments are identified by the argument name in

upper case (not underlined)

.

2-1.22
Section 1 provides a summary of the initial estimates and control values.

It lists the following information.

• The initial values of the parameters, PAR, and whether they are to be

held fixed or not, IFIXED.

• The scale values, SCALE.

• Either the step sizes used to approximate the derivatives numerically,

or, when user-supplied (analytic) derivatives are used, the results of

the checking procedure; and the control values used in these

computations as applicable [see §E.2.b and §E.3.b].

• The number of observations, N.

• The number of observations with nonzero weights, NNZW.

9 The number of independent variables, M.

• The maximum number of iterations allowed, MIT.

« The maximum number of model subroutine calls allowed.

• The two convergence criteria, STOPSS and STOPP.

» The maximum change in the parameters allowed at the first iteration,

DELTA.

fl> The residual sum of squares computed using the starting parameter

values.

» The residual standard deviation computed using the starting parameter

values, RSD.

Section 2, lists selected information about each iteration and includes

the reason the iterations were terminated. The information provided for each

iteration includes the following.

• The iteration number.

9 MODEL CALLS : the total number of times since execution began that the

user's model subroutine has been called, not including calls required

to approximate the derivatives numerically.

9 RSD: the residual standard deviation computed using the parameter

values from the current iteration.

• RSS : the residual sum of squares computed using the parameter values

from the current iteration.

2-1.23
© REL CHNG RSS : the relative change in the residual sum of squares

caused by the current iteration.

• FORECASTED REL CHNG RSS : the forecasted relative change in the

residual sura of squares at the current iteration, and whether this

value was checked against STOPSS ( CHKD = Y) or not ( CHKD = N)

.

• REL CHNG PAR : the maximum scaled relative change in the parameters at

the current iteration, and whether this value was checked against

STOPP ( CHKD = Y) or not ( CHKD = N)

.

• CURRENT PARAMETER VALUES : the estimated parameter values resulting

from the current iteration.

Section 3 provides the following information for each observation,

i = 1, . . .

, N, based on the final solution.

• ROW: the row number of the observations.

• PREDICTOR VALUES : the values for up to the first three columns of the

independent variable matrix, XM , not including the first column if it

is constant.

• DEPENDENT VARIABLE : the values of the dependent variable, Y.

© PREDICTED VALUE : the estimated predicted values, PV, from the fit.

• STD DEV OF PREP VALUE : the standard deviations of the predicted

values, SDPV.

m RESIDUAL : the error estimates, RES.

• STD RES : the standardized residuals, SDRES.

• WEIGHT : the user-supplied weights, WT, printed only when weighted

analysis is performed.

Section 4, displays the following plots of the standardized residuals.

» The standardized residuals versus row numbers.

• The standardized residuals versus predicted values.

• The standardized residuals versus the standardized residuals lagged by

one observation, that is, SDRES(i) versus SDRES(i-l) , i = 2 . .

.

, N.

© The normal probability plot of the standardized residuals.

Section 5, summarizes the following information about the final parameter

estimated and their variances.

2-1.24
• The variance-covariance matrix, VCV, of the estimated (unfixed)

parameters, and the corresponding correlation matrix,

r jk = VCV(j,k)/(vCV(j,j) VCV(k,k)) 1/2

for j = 1, ..., NPARE

and k = 1, . . .

, NPARE.

• PARAMETER : the final value of each parameter, PAR(k), k = 1, ..., NPAR.

• SD OF PAR : the standard deviation of each estimated parameter,

(vCV(k,k)) 1/2 for k = 1, ..., NPAR.

• RATIO : the ratio of each estimated parameter to its standard

deviation,

RATIOk

= PAR(k)/(vCV(k,k)) 1/2 f or k = 1 , ..., NPAR.

• APPROXIMATE 9 5 PERCENT CONFIDENCE LIMITS : the lower and upper 95

percent confidence limits for each parameter, computed using the

appropriate value of the Student's t distribution with NNZW-NPARE

degrees of freedom.

• the residual sum of squares, RSS(3).

• the residual standard deviation at the solution, RSD.

• the residual degrees of freedom, NNZW-NPARE.

• an approximation to the condition number of the derivative matrix, D

(the Jacobian) , under the assumption that the absolute error in each

column of D is roughly equal. The approximation will be meaningless

if this assumption is not valid; otherwise it will usually

underestimate the actual condition number by a factor of from 2 to 10

[see Dongarra et al. , 1979, p. 9.5].

NOTE : The standard deviation of the predicted values, the standardized

residuals, the variance-covariance matrix, the standard deviations of the

parameters, and the 95 percent confidence limits on the parameters are all

based on a linear approximation to the model in a neighborhood of the

solution; the validity of this approximation depends on the nonlinearity of

the model. The statistics based on this approximation may be extremely

inaccurate for a problem with a highly nonlinear model.

E .2 Step Size Selection

E.2.a The Step Size Selection Algorithm

The STARPAC step size selection subroutines use an algorithm developed by

Schnabel [1982] to compute optimum step sizes for approximating the partial

derivatives of the model with respect to each parameter. Briefly, the

2-1.25
relative step sizes selected by these subroutines are those which produce

forward difference quotient approximations to the derivative, Df^j that agree

reasonably well with the central difference quotient approximations, DC(j. The

central difference quotient approximations are twice as accurate but also

twice as expensive to compute. Since the additional accuracy is not usually

needed, central difference quotient approximations are not used by the

estimation subroutines.

The number of reliable digits in these derivatives is a function of the

step sizes used to compute them. Given properly chosen step sizes, the number

of reliable digits in Df^ and Dccj will be approximately n/2 and r\ , respectively, where n is the number of reliable digits in the predicted

values, PV, from the user's model subroutine. For example, if the predicted

values are computed using an iterative procedure (such as quadrature or a solution of partial differential equations) which is expected to provide five

good digits, then n would be five; if the predicted values are calculated from

a simple algebraic expression translated directly into Fortran code, then n would (usually) be the number of decimal digits carried by the user's computer

for the results.

The relative step size for B(k), k = 1, ..., NPAR, is initially

STP(k) = 2(1CTNETA

/ Y ) 1/2

for k = 1, ..., NPAR,

where

Y is the average curvature (estimated by STARPAC) of the model with

respect to 3(k)

.

NETA is the argument used to specify the number of reliable digits, n , in

the model results.

The forward difference quotient approximations with respect to 3(k),

k = 1 , . .

. , NPAR are then

f

i (x i ,Bk

) - £± (x±i $)

D fd (i,k) = for 1-1, ..., N,

STP(k) -SCALE(k) •SIGN( B(k)

)

where

N is the number of observations.

f^ is the function which models the i tn observation.

Xj is the vector of the values of the M independent variables at the i

observation.

B is the vector of the NPAR parameter values.

th

2-1.26
8 is a vector which has the same values as 8, except that the

k

1

- parameter is equal to

8(k) + (STP(k).SCALE(k).SIGN(3(k))).

SCALE is the vector of the scale value (i.e., of the typical size) of each

of the parameters.

SIGN is a function which returns the sign of its argument.

The central difference approximations to the model derivative with

respect to B(k), k = 1, ..., NPAR, are

f

i (x ±i B+h - f± (xlf B-k

) Dcd (i,k) = for 1 = 1, ..., N,

3 1/3 .10"NETA/3 .SCALE(k).SIGN(B(k))

where

8 is a vector which has the same values as 8, except that the k*- parameter is equal to

8(k) + (3 1/3 .10"NETA/3 .SCALE(k).SIGN(8(k))).

—k t" h

8 is a vector which has the same values as 8, except that the k

LU

parameter is equal to

B(k) - (3 1/3 .10"NETA/3 .SCALE(k).SIGN(8(k))).

The relative step size is considered acceptable if, for at least N-a

observations,

|D fd (i,k) - Dcd (i,k)| < min{lO-NETA/4 )10-2| for i = 1, ..... N,

where a is the number of observations exempted from meeting the above

acceptance criterion [see §D, argument EXMPT] . If the step size is not

acceptable, it is adjusted by factors of 10 until the condition is met, or until no further decrease in the number of failures can be made, although in

no case will the selected relative step size be greater than 1.0 or less than

j^-NETA^

Note that the step size selection subroutines will return the selected

step sizes even when the number of failures exceeds the allowed value; this

condition will be noted by the value of IERR. The detailed printed output

should always be examined for problems discovered by the step size selection

subroutines.

2-1.27
E.2.b Computed Results and Printed Output

From the Derivative Step Size Selection Subroutines

The argument controlling the printed output, NPRT, is discussed in §D.

The output from the step size selection subroutines consists of a summary

of the input and control values and, for each parameter, the selected relative

step size, the number of observations at which this step size failed the step

size selection criteria, and the row numbers at which the failures occurred.

E .3 Derivative Checking

E.3.a The Derivative Checking Algorithm

The STARPAC derivative checking subroutines use an algorithm developed by

Schnabel [1982] to determine the validity of the user-supplied derivative

subroutine. The user-supplied derivative subroutine is considered correct for

a given row i, i = 1, . .

.

, N, and coefficient B(k) , k = 1, ..., NPAR, if

D fd (i,k) - D(i,k)| < 10-t |D(i,k)|

where

'fd

is the derivative computed by the user's subroutine.

is the forward difference quotient approximation to the derivative.

[See §E.2.a.]

t is the agreement tolerance, i.e., number of digits of agreement

required between D and Df^, which must be less than or equal to the

number of good digits in Df^. [See §D, argument NTAU. ] When the agreement tolerance is not satisfied, the checking subroutine

attempts to determine whether the disagreement is due to an error in the

user's code or is due to the inaccuracy of the difference quotient

approximation, caused either by high curvature in the user's model or by

significant roundoff error.

The derivative checking subroutines each check only one row of the

derivative matrix. The user should examine the row at which the derivatives

were checked to ensure that some relation between the coefficients and

independent variables, such as a zero-valued independent variable or a factor

(x^ - 3(k)) when x^ = $(k), is not hiding the effect of an incorrectly

computed derivative. Checking only one row is appropriate since the same code

is frequently used to compute the model function and derivatives at each row

i =1, . .

. , N, as is the case in the examples shown in §F. If the code used

to express the model function and derivatives is not the same for each row,

then each distinct section of the code should be checked by making multiple

calls to DCKLSC with NROW set to a row within each section. [See §D, argument

NROW.]

2-1.28
E.3.b Computed Results and Printed Output

From the Derivative Checking Subroutines

The argument controlling the printed output, NPRT, is discussed in §D.

The output for the derivative checking subroutines consists of a summary

of the input and control values, and the results of the derivative checking

test with respect to each of the model parameters, 3(k), k = 1, •••, NPAR.

The possible test results are:

OK -

• The user-supplied derivative and the numerical derivative agree to the

required number of digits.

QUESTIONABLE -

• The user-supplied derivative and the approximated derivative agree to

the required number of digits but both are equal to zero. The user

should recheck the derivative at another row.

• The user-supplied derivative and the approximated derivative do not

agree to the required number of digits, but the user-supplied

derivative is identically zero and the approximated derivative is

approximately zero. The user should recheck the derivative at another

row

• The user-supplied derivative and the approximated derivative disagree,

but the user-supplied derivative is identically zero. The user should

recheck the derivative at another row.

o The user-supplied derivative and the approximated derivative disagree,

but the validity of the approximated derivative is questionable

because either the ratio of the relative curvature of the model to the

slope of the model is too high, or SCALE(k) is wrong.

• The user-supplied derivative and the approximated derivative disagree,

but the validity of the estimated derivative is questionable because

the ratio of the relative curvature of the model to the slope of the

model is too high.

INCORRECT -

• The user-supplied derivative and the approximated derivative disagree,

and there is no reason to question the accuracy of the approximated

derivative.

2-1.29
F . Examples

The sample programs of this section use the model and data given in

example one, pages 428 to 441 of Daniel and Wood [1980]; the model is

f± (x lt 3) = 6(l)-x1 (l)3( 2

), for i = 1, ..., N.

In the sample program of figure F-la, NLS is used to compute the least

squares solution using numerically approximated derivatives; figures F-lb

through F-lf show the output from NLS.

In sample program of figures F-2a, NLSD is used to compute the least

squares solution given analytic derivatives; figures F-2b through F-2f show

the output from NLSD.

In the example program of figure F-3a, STPLS is used to compute the

optimum step sizes for numerically approximating the derivatives with respect

to each of the parameters, B(k) , k = 1, 2. Figure F-3b shows the output from

STPLS.

In the example program of figure F-4a, DCKLS is used to check the

validity of a user-supplied derivative subroutine. In figure F-4a, the

derivative subroutine has been intentionally coded incorrectly in order to

display the report obtained when the derivative checking subroutine determines

the derivatives are incorrect, and the starting parameter values have been

chosen in order to display the report obtained when the test results are questionable. Figure F-4b shows the output from DCKLS.

2-1.30
MAIN PROGRAM!

MODEL SUBROUTINE!

DATAl

Nonlinear Least Squares Estimation

With Numerically Approximated Derivatives

PROGRAM EXAMPL (INPUT, OUTPUT» TAPES INPUT. TAPE6 OUTPUT)

REAL YUO), XM(10»5), PAR«5), RES(IO) DOUBLE PRECISION DSTAK(2O0) EXTERNAL NLSMDL COMMON /CSTAK/ DSTAK OEFINE VARIABLES FOR VARIOUS DIMENSIONS

LOSTAK 200

IXM 10 READ NUMBER OF PARAMETERS STARTING VALUES FOR PARAMETERS NUMBER OF OBSERVATIONS AND NUMBER OF INDEPENDENT VARIABLES INDEPENDENT AND DEPENDENT VARIABLES

READ (5, 100) NPAR

REAO (5, 101) (PARII), I»1,NPAR)

REAO (5, 100) N, M READ (5, 101) MXH(I,J), I-1.N)# J"1,M), (Yd), I-1,N) PRINT TITLE AND CALL NLS TO PERFORM NONLINEAR REGRESSION

WITH NUMERICALLY APPROXIMATED DERIVATIVES

WRITE (6, 102) CALL NLS (Y, XM, N, H, IXM, NLSMDL, PAR, NPAR, RES, LDSTAK)

STOP FORMAT STATEMENTS

100 FORMAT (215) 101 FORMAT (6F6.3)

102 FORMAT (1H1, 44HTHIS IS AN EXAMPLE OF STARPAC SUBROUTINE NLS/)

END SUBROUTINE NLSMDL (PAR, NPAR, XH, N, M, IXM, PV) SUBROUTINE TO COMPUTE PREDICTED VALUES OF OEPENOENT VARIABLE

REAL PAR(NPAR), XH(IXH,M), PV(N)

DO 10 1-1,

N

PV(I) - PAR(l) * XM(I, 1) • PAR(2)

10 CONTINUE

RETURN

END

0.725 4.000

6 1 1.309 1.471 1.490 1.565 1.611 1.660

2.138 3.421 3.597 4.340 4.662 5.660

Figure F-la

Example program and data using NLS

2-1.31
k-i UJU

<7 I—

u- o u. » UJ • » •» * » "- » » I— » * * • * » * • — • • at * » uj » • o • * * • o * • UJ • » >- » » * * » JT * * >- » • x * » O • * a; * » a. • # o. * * < * C c— I— u">

1/1w • « •» » 1_J • I— * • a. • • uj * » x • » 3 * » z * » * • X • » »- » • —• # » -» » • • • Z » • o • c QC -J UJ

QCO *T c -J *- o z UJ O o o e—

U- l~< :> 1 1 u. QC <W M UJ UJ O uj <J t— ^ <\J f» UJ

isi

r •< a. CD (^ I Z

•—* « > t— o f- fC

1/1M t— <•> 00

m ^-»

a at

w

m •-« t/t Ja- Q. « UJ o er* < il

uj a.

a 0*

O I UJ

i— a. »f*

m—

I

0)3

t/>

t-1 l/l 3

UJ —J I «» 3a

•u3Owz

CO w

-i-i 3 3 o «-> u. u.

|/| (/^ UJ UJ — a a • UJ * z * » a * 1/1 • b— • JJ • - * at * •—

» « * a • 3 » z » O * o

1/1 » o < » -I • x » •

UJ • QC • » •« * Z • UJ •

CO * —i »

«- • z • * o *

1/1 * Z • z — o o — a. o O

*- * o o a a o o « — in o - i\j o

i/i r- •

UJ

C

X UJ

O «/» I t_ 3T •— z - • t/i UJ

t- o Q V\ z X « •- z UJ t- o UJ > - I-* t— ImU •—

a. a. o o - 1^ o^

UJ JX UJ z e>-e 4 z z

VI UJ oc

O » o o

CS

M 3 t/i QC »-d I—

O UJ

•< QC —J «/> « •«

U. u. > 3 a » » a o QC ^ << o QC oc 3 ** •—

UJ UJ

a: oc

u « _» u> v\ i/i UJ UJ u. UJ o CD to

<r> a>T QCO o z Z O -J z 3 3 —

U, o u_ u. z z X « o

1-

OO 1 1 at QC at at UJ o UJ UJ o z « 1^ z ac 1—

—i -i CD <f < 3

—1 b- UJ

X a.

a

-« 3 »— 15 o Z o (/I Z Ul

oo at UJ < uj

X z to H- X > - at a 3 o »-* o

** v> oc >- z u. C U.

M

aa * 3 •4

X. i-i -1

a o 3

Of UJ UJ -J o

u. UJ oc -J —J tSi O f- « •< ** o c> U. at at UJ 1/1 UJ o

UJ o H- o

CD t/> z z z

3. UJ 4 3 < 3 3 QC 9 UJ 3 t

2-1.32
ac

a >- •a x

a.

T

<_>

O

H

sr o T 1 U uj ui n o

_j _j rr UJ <3 fQt > -^

vo> « sr a. T

«_> IS r- y O T I u UJ uj

_i _i e

OI » i-<

0)3

</>

O

V

C

w

X

uja 1

»- o • T

<-j I UUUJ O -I _i o C

w

X

uja

X - o

ij*) tj5 l\J « z »< o I I UI U UJUJ C -j —I -o U. UJ < f*- ar =» JD 1 3

Ci- &

4-1 aj 3

)-i

O

a C/i

•i-i ,J ^ Zae

(0xw

UJ =>

IO _l

ir z «

OJ .- >

&XUJ UJ - O -J UJ

_J V* CTs

O _) r o

2-1.33
coat •••••• H M ^ r^ ^s -5 o o o o o 1 I I I I I tXJ m '" Ul UiUI

_J f\|

O <—• f\J •-* <NJ <3 {> MO^OO 3 ^ GD <-* ® **- -f a n o a o ^ o

w=i pH rf*» *X>KT CI

i/>

—# (O ir> -c aD UJ ^ <B (\l t+* -O <0 at m O h r* rflfl I I C 3 -X <E (NJ .* »H <M

_jO lO ^ OD HK

:> « o u ifi^ Mm

UJ "> P*1 -O r-#

O ** flO

OO -*

OO IT rt

O rvj £ ^ 'T <£> 'O

O lii f\j •* ** ^ ^* ^*

t— o£ •••••• c o *r\ ^ <r ^ *^ "*

uj r- -r

O *^ f* *o

h* uj

h ^ ^*' Ort

(J 3 rl

H * ^3

O *

1-..J ^T H-# (\J IT 4)

Qt ••••••

a. (\j in

m j^ # ^

4-1 coV3

t— uj OOOOOO Z-J

ujcB OOOOOO

o< OOOOOO 3. _-» OD rfNOMO

ujQt

w r\j^ * (D-O

Q. •<

H -t IT <*! ffl

O

O

M fO (« * * if>

0) 3O&e

OOOOOO

OOOOOO o *< o ir w o o f- a -o t-< «

en -r *» tTt .0 -o

2-1.34
CD I I I I » I I 1 1 * I I I I <• I I < i * I J I ! * ; i I i * I i I i * i I i i i i I i * i I i I * ^

O UJ

*/l O

a

llll*llll + lltl*!lll+l!ll+<-» Itll+llll + llll+llll + llll

III lf|lll*lllt+lltl*llll*0 + I I I I * I I I I I I I I + I I I I * i] 1 I I

I 0)3

be

•H

« or *

»- (y.

QC CO

0)3c

3

a

4-1 3

o

,J2

+ llll*llll*llll+llil»ill r- t|ll*llll*llll*llll + lltl+* «% if\ m M\ M> m «w

t«. «m r- ^ «m r» i

2-1.35
r a. -l «si —

a. *\J

O

_j => o o

UJ -> o o <s •

UJ z

UJ • » < Q r ~ d

•HUa

l— ct

a -o Z lu mo

uj 1 m©

t_> o o rat -J <r-» -o UJ <e c-» 1 3

fa

a

4-1 0) 3

J-l

O 3M CO

•H kJ

fa Z o yO 1*) (VI ^

»-• o -o o o

>— • • 1 1 •« CVJ *• UJ UJ

at. * r» a> vH o r*

aE

cC

fa

X ^ («1 CSJ ac f^

r

t— o o 4 o o

_l 31 1 1 a. 0»

n 2 « o UJ UJ w* o 1 o z _J - »>» u. gg CM — o UJ •o m o f\l l>- O

»- o GO f -e CD—

a. < —

o> ir>

a f*

m 1 >- •« UJ IM o r- en • • ~ Q z ac m vo I—

ts\

a « 0> f\J

M

ist UJ o U_ « X < CO 1 - t—

(- v> ct z o X LU z z Ui C UJ UJ at o a cr > UJ »-» « «- Q X. c X <-> CO 3 or J3 o I- UJ 3

UJ cc >- t-m o O UJ CM in </« « UJ z

V. « U. • uo t— <M o UJ

— ac < r U. UJ UJ •o •r ac > u. z

CO UJa UJ l*> 0> »—X GO o « UJ

a

ata o ^ en 4 CD •o 3 Q U. »-• z •* UJ u tn r» 4 at •O CO a a - a ac »\j r. UJ 4 r» • t/> a -< •* VI <i z «-4 •*• O —J a. • m at i/> o

1- UJ o m o> u. < UJ z

•« «_l ts\ *- «n 0» z o O UJ o Z Z UJ - • • a z ac o — « <-) a 1 IX a z <« o X •— z —i u- UJ c? a 3 t- UJ UI O Or. <* UI X z z VI V> Q (- a < »—

ox i/» •-• «

a. > ac IX UJ u. «J _/ z JT a O « o 1— 4 <« a •-• •« <_> =» <_> z J* «M •* 3 3 >» r z *t O o a o • 1 1 1 3 •-* UJ •^ •« u* ac a * «u

2-1.36
Nonlinear Least Squares Estimation

With User-Supplied Analytic Derivatives

MAIN PR06RAMJ

MODEL SUBROUTINEi

DERIVATIVE SUBROUTINE

OATAl

PR06RAM EXANPL (INPUT, OUTPUT, TAPES • INPUT, TAPE6 • OUTPUT)

REAL Y(10)» XM(10,9), PAR(5), REStlO) DOUBLE PRECISION DSTAM200)

EXTERNAL NLSMDL, NLSDRV COMMON /CSTAK/ DSTAK DEFINE VARIABLES FOR VARIOUS DIMENSIONS

LDSTAK • 200

IXM - 10 READ NUMBER OF PARAMETERS STARTING VALUES FOR PARAMETERS NUMBER OF OBSERVATIONS AND NUMBER OF INOEPENOENT VARIABLES INDEPENDENT AND DEPENDENT VARIABLES

READ (9, 100) NPAR

READ (9, 101) (PAR(I), I-1,NPAR)

READ (9, 100) N, M READ (5, 101) ((XM(I,J), I-1.N), J-1,M), (Yd), I-1,N) PRINT TITLE ANO CALL NLSD TO PERFORM NONLINEAR REGRESSION

WITH USER-SUPPLIED DERIVATIVES

WRITE (6, 102) CALL NLSD <Y, XM, N, M, IXM, NLSMOL, NLSDRV, PAR, NPAR, RES,

* LDSTAK)

STOP FORMAT STATEMENTS

100 FORMAT (215) 101 FORMAT (6F6.3)

102 FORMAT (1H1, 44HTHIS IS AN EXAMPLE OF STARPAC SUBROUTINE NLSD/)

END SUBROUTINE NLSMDL (PAR, NPAR, XM, N, M, IXM, PV

)

SUBROUTINE TO COMPUTE PREDICTED VALUES OF DEPENDENT VARIABLE

REAL PAR(NPAR), XM(IXM,M), PV(N)

DO 10 I"1,N

PV(I> • PAR(l) XM(I, 1) • PAR(2)

10 CONTINUE RETURN

END SUBROUTINE NLSDRV (PAR, NPAR, XM, N, M, IXM, D) SUBROUTINE TO COMPUTE THE PARTIAL DERIVATIVE (JACOBIAN) MATRIX

REAL PAR(NPAR), XM(IXM,M), 0(N,NPAR)

DO 10 I-1,N

0(1,1) - XM(I,1) ** PAR(2)

D(I,2) PAR(l) * XM(I,1) ** PAR(2) * AL0G(XM(I,1) ) 10 CONTINUE

RETURN

END

0.729 4.000

6 1 1.309 1.471 1.490 1.969 1.611 1.680

2.136 3.421 3.997 4.340 4.682 9.660

Figure F-2a

Example program and data using NLSD

2-1.37
o <r o O -H « o »- * < * * > * * •— » * rx » * LLi * * e ** o * » uj » » ^ » * _l • * Q. * * a * * 3 »

42CNI

QJ

l-i 3

•H 3&

•u3Q

•J

1Z

14-) o

CO w * =3 * o * # T * * y— * * t— * * 3 * o -3- * 7 * V a -& » <—

* » h- * * <1 * • r * * t—

* * fr- » « iz' ft * UJ * * & * i^-i » • U_. * * a * * -* » * _ * • o * * v * * * & (— O * ^ * « « * • Li-. * « 2 » * O » • z •

u.

a

or uj I «- - - «> o *

—i —t Z> 3 _) UJ ' a « •« « 3 15 X O U U. u. 1/1 Z O I !/)</) UJ UJ UJ »—

I *"- O O CC 3£ UJ . u or i _j UJ UJ ' UJ X Si e o i © © - » —• : O O «— a at i © © O Uj U. I O © •-• O O I o © o l ir o z x : rj © uj »-» 5ji

^- • -J *-4 • * CD Irt XI C C z z x UJ

UJ 3 O -J c </» < H- > i-

_J

a LU

—1 UJ _j - or

T « X < UJ < UJ «

t- 3 X »- Qt

</1 o - ts\ UJ <a —J 1 t-4 a. X a

-J c •/> z UJ «

«3 UJ J— t- OC -

<_J o Of UJ < => UJ UJ z. a o. UJ */1 z O * z z « t-« z IX >- !-•

l—l a> < 4 3

»- UJ

X a.a

a 3 i— 19 CJ z a o i/1 z UJ t—

u.

or UJ < UJ

X

<r tX > i— or z

r> *_> •-H

n o

4J1 tt. t- z u. !* a UJ < N* (- -J u. > ^ i/> «

UJ !-• UJ o UJ 1—

o z ^ IX UJ or > o a « 3* « UJ 4 3 -

2-1.38
UI1IUJ 3 O

_| -J »

UUiUJ

_l

—I .O uj <t O

T3

0)C

</i

av

O l/lX

oj tx I

(UU UJ

CX 3 1* Q -l _l O

ISI o >-

O </>

X

UJ OL I - «_)

i/» CD

N « z i-t O I 1 UUtb UJ oI

0)

s-l 3a

4-1 3Oo

Ujz

1/1

M

t/> ** a 1 uu

(J r»- Z o I v?

uj 3

IQ _J z z •*

C

X UJ UJ 3

tO _J

2-1.39
o ir ^ m o o

•t er% *r

m i\j co

<_i <\j

—« <v *< *-* o o o o o o I I I I I I lil Hi uj Hi UJ 111 _l o **• <\' **• o «* < 0> <* iT« m O 0* 3 * oo -* co f^ wr O N OP O M «0

«—j H IT 3} £ i>

n

</* r-i ^ irvm -O <D LU ^ CO <\J ff* •/* «0

at (*i O* ^ ^» <*> l*i I I o o o o o 1 I I I I I \±, l±i ^J ^ ^J I^J LLUJ *- lf\ i-l *" <NJ f-

_l Of f O^^1 UJ> f- <T> —i £> — CD O O'T'CO^ h o r\j o ri <r /) o

quj rvi **-*#* ^ *sj t— a ••••• i/i a o m o> so & m »r> lu ^ «#- o ** *- -O

i— uu ^ ir fj * o w u D h h^ £m» <-t j .* ^ -J- «\i ir ^>

QC ••••••

a. cvj rnm <r *# *»

*- lu oooooo

Z_i oooooo lu CC OOOOOO O < OOOOOO Z •— CD ^i ^ O M O

lu ctm MP * » *

T3I in

0)

4-1 3a

4-1 3OcJ

M—

o

ro I «« » ITi

cC w

_i c c o o o o

•a o o o e O o » oooooo OOOOOO O O f* 9 -C »H <x> H- m -r -e *> <c <o

2-1.40
Illl*llll*llll*llll*lllt lllt*llll*llll*llll*!lll*«>

at*

c: <u3

C

•III

l * l i i i l l

i l i

i i i i i + i i i

l **

I o

•

•a

i i i i iii*iii i i i i * i i i i *

*l|ll*llll*ltll*llll*llll*m

CD 4J CN 3

1 a

t* 4J3

01 o

u

3 Q

bC c/> •H J

b z

U-lo

£

CO w

"" or •

*iiii*iiii*iiii*iiii*iii

o

•

i ** i l i l iiii*iiii*iiii*iiii*« % ir* ift *» ** «*»

2-1.41
0)3C

>« z C D a <-> a.O 1- ffi « IT 4ZIU m o» uj 3 <*> ao

<_>

C

© fa —i

«4 <C LU CD e*

CM s 1 a

Eu 4-1 3

0)

O ^s a

w: t/>

O *> en

1-4 o >o oI o1 LU UJ

ec f-i <#• O **

el

m

f- in

f-a CO m <M •

.Jo&

z. < o O Z -I

i-t

O 11/ h- o cr

et « _i

T" •— -*3 LU o C «

if u. o < T < *^

H- -• - Z O Z m (V O O

u< LU o <\J f^ •# Kl -o o «r o f~

<•» <c

cy <M o c T O m Z>

1U CD )— •— o o VI * U. 1 </l « Z u- UJ

CD LU

C UJ c o Hat a O •* V) z •« UJ o fr> ^ <*J c Or «m r» UJ

l-< l/» « z t-# » o -J

*- LU o m o> « O </l —• m o> z z z LU t— • • o c— « <_> « • ec

•m i-* z -J

u

Ca <a UJ

or « t—i a. v>

a. > OL c LU

OC « a »- ) 1 >t I zs3

-Jo

—1 <M

r

Ot -T f-< * ^ fi

a. c .o — -O

U. CO fvj C M l>- CO _j o ** in

a c CO

—9 c JE o © - UJ r> ej m o «I UJ z

«M

O Ui *-<a c * ot > u. z

co o « LU c

CO ^ =>O LL »-» •C CO o O tr- • </> a -• • tn a «/» o

u. « Ui z o o Ui o z Ot o z •« iS

OO 3 (- UJ UJ z z tn n o

_J—1 z I « •a o ^» 3 a K O o o o

»-t *-« UJ

a

fi rvi «/•> on « o. UI UJ o a. at at a 4

2-1.42
Derivative Step Size Selection

MAIN PROGRAM PROGRAM EXANPL (INPUT, OUTPUT* TAPES • INPUT, TAPE6 • OUTPUT!

REAL *H!lO»5l, PAB(5)* STPI3) DOUBLE PRECISION DSTAM200)

EXTERNAL NLSNDL COMMON /CSTAK/ DSTAK

C

C DEFINE VARIABLES FOR VARIOUS DIMENSIONS

C LDSTAK - 200

IXM 10

CC READ NUMBER OF PARAMETERS

C STARTING VALUES FOR PARAMETERS

C NUMBER OF OBSERVATIONS AND NUMBER OF INDEPENDENT VARIABLES

C INDEPENDENT VARIABLES

C READ (5, 100) NPAR

READ (5, 101) (PAR(I), 1*1, NPAR) READ (3, 100) Ni M READ (3, 101) ((XM(I,J>, I-l,N), J-1,H)

C

C PRINT TITLE AND CALL STPLS TO SELECT STEP SIZES FOR

C APPROXIMATING DERIVATIVES

C WRITE (6, 102) CALL STPLS (XN, N, H» IXM, NLSMDL, PAR, NPAR, LDSTAK, STP)

C STOP

CC FORMAT STATEMENTS

C 100 FORMAT (215)

101 FORMAT (6F6.3)

102 FORMAT (1H1, 46HTHIS IS AN EXAMPLE OF STARPAC SUBROUTINE STPLS/)

END MODEL SUBROUTINE! SUBROUTINE NLSMDL (PAR, NPAR, XM, N, M, IXM, PV)

C

C SUBROUTINE TO COMPUTE PREDICTED VALUES OF DEPENDENT VARIABLE

C REAL PAR(NPAR), XM(IXN,M), PV(N)

C DO 10 I 1, N PV(I) • PAR(l) • XM(I, 1) •* PAR(2)

10 CONTINUE

C RETURN

END

DATAl 2 0.725 4.000

6 1 1.309 1.471 1.490 1.565 1.611 1.680

Figure F-3a

Example program and data using STPLS

2-1.43
3 o ct -J U.' u.

at is -#• -o .J Ha c z UJ © o a i—

IX 1— 5» • 1 u_ Ot UJ - t-« UJ Uj O to «

UJ ^ *— ** •o <? UJ 3 t-» w z. <3 a. o •> I z « at f* !— » t- O IT 1- o u UJ

to J* t~t to <\J r» t—» 1— « O C£ ** •O OD to *- UJ »* »-» a. or UJ O O « t_> X at a

uj a. o o r- I UJ l— (_i UJ

>-

a <C w>

ato

—Jto

UJ

INI o

UJ zc

1—

t-Z

u. •-• f— UJ o

COI

a) 3O

-J -I 3 3 « « « « O O U u. to to UJ u,' • r- » • D * * *« • » - * UJ » <_> » 3

UJ * -J * _J • a « • UJ • UJ > » to • H- o o • • UJ IS «-. o o * u. • t 2 a. o o • fj • « >- « o o • ft • IX •- o. o o » to » « Qt *^ IT C » • O. « <M O • a. • - r- • * UJ • to • •# * C£ • * UJ * * O •

to o —J •-• UJ UJ t—

a- a. to to tj

UJ to •— UJ

t- 3 -j

x ->

to to 3 O UJ

to at to

19 to UJ U. Z • ** a

X < o e C —J UJ —

—J UJ Of »— ji uj uj ;— u- •< O O DO. u I O X O

_i

x z uj uj

to s x t- z uj z uj a. a to x •— z

*-i

T *— • to UJ

»- o o to z x

•« - Z UJ »-

O uj

j> t- Ml- t-l

K

at o. o tJ> •— to to ur uj z •-' * z z

to uj at o > o o

cdx 3 to

a «-• »— PU »- Uj uj uj »— •— U.U. >3 CD C > > OO C£ —I « O Ot Ot S M M UJ UJ

OC Of t_>4 _J u. to to

uj uj u. u* o m er

a-> cc

X or

OO Z X O -J Z 33 *- —j u. a u- u. zz x « o •-> a c

*- I I at ac at ae UJ O UJ UJ

u. o to O.

W CO

X

OX

X 3 at 3 3 z »• z x

2-1.44
Derivative Checking

MAIN PROGRAM!

MODEL SUBROUTINE!

DERIVATIVE SUBROUTINE!

PROGRAM EXAMPL (INPUT* OUTPUT, TAPES INPUT, TAPE6 - OUTPUT)

REAL XM(10,5), PAR(5) DOUBLE PRECISION DSTAK(ZOO) EXTERNAL NLSMOL, NLSORV COMMON /CSTAK/ DSTAK DEFINE VARIABLES FOR VARIOUS DIMENSIONS

LDSTAK 200

IXM « 10 READ NUMBER OF PARAMETERS STARTINS VALUES FOR PARAMETERS

NUMBER OF OBSERVATIONS ANO NUMBER OF INDEPENDENT VARIABLES INDEPENDENT VARIABLES

READ (5, 100) NPAR READ (3, 101) (PAR(I), I- 1, NPAR) READ (5, 100) N, M READ (5, 101) ((XM(I,J), I-1,N), j>l,M) PRINT TITLE AND CALL DCKLS TO PERFORM DERIVATIVE CHECKING

WRITE (6, 102) CALL DCKLS (XM, N, M, IXM, NLSMDL, NLSORV, PAR, NPAR, LDSTAK)

STOP FORMAT STATEMENTS

100 FORMAT (215)

101 FORMAT (6F6.3)

102 FORMAT (1H1, 46HTHIS IS AN EXAMPLE OF STARPAC SUBROUTINE DCKLS/)

END SUBROUTINE NLSMDL (PAR, NPAR, XM, N, M, IXM, »V) SUBROUTINE TO COMPUTE PREDICTED VALUES OF DEPENDENT VARIABLE

REAL PAR(NPAR), XM(IXM,M), PV(N)

DO 10 I-1,N

PV(I) " PAR(l) XM(I,1) PAR(2)

10 CONTINUE RETURN

END SUBROUTINE NLSORV (PAR, NPAR, XM, N, M, IXM, 0) SUBROUTINE TO COMPUTE THE PARTIAL DERIVATIVE (JACOBIAN) MATRIX DERIVATIVE WITH RESPECT TO FIRST PARAMETER HAS BEEN CODED INCORRECTLY (WITH MULTIPLICATION RATHER THAN EXPONENTIATION)

FOR DEMONSTRATION PURPOSES

REAL PAR(NPAR), XM(IXM,M), DIN, NPAR)

DO 10 I-1,N

D(I,1) XM(I,1) • PAR(2)

0(1,2) • PAR(l) » XH(I,1) * PAR(l) * ALOG(XN(I,D)

10 CONTINUE

RETURN

END

DATA! 0.000 4.000

6 1 1.309 1.471 1.490 1.365 1.611 1.680

Figure F-4a

Example program and data using DCKLS

2-1.45
« 3 3

t- « O

Ul - <x Iu

P4J3O

O. Z

QC

O

a «

<_> »- z </->

PK

CO

— »- »- -i •<

a

-J _i 3 3 O i-

*i * « «i u.

a a:

O

O U. U. UJ Ui

V> l/l UJ Ui Ui OX — O Q I tUi «I Xu «/» « Ui u. a x z tart i-i -

a. 3 Ui at o tart *- tX

at -I at a. ui •* rtf Ui « (X « a. X </> > » o o - Ui > o to « o - tart " z o

wi * » » H-

— o Ui Ui tart at OC Ui o • * WlflK o </» O OS u Ui Ui ft. & «/> u. » o » r z •* o Ui z tart a O UI © z a » z » < « a. o X •I a e m o * -• * KK~ o *- • z XX** • tart Ui • * * 4 at o z o o

-j » <_> * a. «i • Ui lu at -j •-* rtj I • X »

tart tart taX Ui Ui

>-! «I u.

O -I « 1 or o o z ««

at -• > o

i/i tt

X u. u- Ui </* u.

a Ul t—

O

a as u>

a

Ui c«

O Z 3 m 3 oo at QC 3 -i at

X z Ul Ui Z * Ui

Ui 3 rttaaa > «o

a HN Z tart

X z » 1 z z ** 3 3 o 3

2-1.46
Appendix A

WEIGHTED LEAST SQUARES

Weighted least squares can be used to eliminate observations from the

analysis, and to compensate for unequal variances in the observational

errors.

Observations can be eliminated from the analysis by using weight values

consisting only of zeros and ones. This will produce the same results as performing an unweighted analysis with the zero-weighted values removed except

that the predicted values, the standard deviations of the predicted values,

and the residuals of the zero-weighted data are computed. There are two main

reasons for weighting observations zero. The first is to obtain the predicted

values and their standard deviations for a set of independent variables not

included in the observed data. (This is done by assigning any arbitrary value

to the dependent variable of the desired set of independent variables, and

then weighting these values zero.) The second reason is to allow easy

examination of the effect of outliers and influential data points. Outliers

often appear as large values in residual plots. Careful checking of the data

often leads to confirmation that the data are in error, and sometimes to a correction. When a cause for suspicious data cannot be found, it may be

advisable to compare the analysis with and without the questionable data.

Caution is in order if the estimates or conclusions are highly sensitive to a small amount of suspicious data. Data that have a very high influence on a fitted curve may not result in large residuals, however, even if they are in

error. In fact, extremely influential observations may force the fitted curve

to be very close, leading to very small residuals. It is therefore desirable

to identify influential observations and to compare the results obtained with

and without these points. Several methods for detecting influential

observations are discussed in Bement and Williams [1969], Cook [1977], Hoaglin

and Welsch [1978], and Belsley et al. [1980].

Using weights to compensate for unequal observational error variances is

not as straightforward as using zero weights to eliminate observations from

the analysis. When the variances of the observational errors, e^, are not

equal, the unweighted least squares estimates remain unbiased but do not have

minimum variance. Minimum variance estimates are obtained by using weights

wt.j = 1/Variance[ e^] when the error variances are known . If weights must be

estimated, they should be based on at least 10 degrees of freedom [see Bement

and Williams, 1969]. In practice, however, weights are derived from theory, or obtained from the data being fit, and either of these methods can do more

harm than good. When the need for weights is suspected and the error variances are not known, first fit the data using unweighted least squares;

analysis of the residuals may confirm the need for weighting and may also

provide estimates for the weights themselves. If the need for weights is

confirmed, then a statistician should be consulted to assist in selecting the

weights and in interpreting the results.

2-A.l
Appendix B

ESTIMATING THE NUMBER OF RELIABLE DIGITS

IN THE RESULTS OF A FUNCTION

The number of reliable digits, n, in the results of a real valued

function, g((3), can be estimated in most cases by evaluating

n = -log10 ( max (lg(B j

) - [a + j»b] I}]

J—

2 fiOOl

where

6 J is the vector of the NPAR parameters of the function given by,

BJ(k) = B(k) + j.lO-( DIGITS/2 ).B(k) for j = -2, ..., 2,

and k = 1, . . ., NPAR,

where

DIGITS is the number of decimal digits carried by the user's computer

for a single precision value when the single precision version of

STARPAC is being used and is the number carried for a double precision

value otherwise.

2 a = (0.20). y g(3J).

J—

2 b = (o.io). y j. g (0j).

j—

This procedure may underestimate the number of reliable digits if g(3) is

extremely nonlinear. A more elaborate and more robust procedure is described

in Gill et al. [1981] .

2-B.l
References

American National Standards Institute. (1966). ANS FORTRAN X3. 9-1966. New

York: American National Standards Institute.

Bard, Y. (1974). Nonlinear parameter estimation. New York: Academic

Press.

Belsley, D.A.; Kuh, E.; Welsch, R. E. (1980). Regression diagnostics. New

York: John Wiley and Sons.

Bement, T. R.; Williams, J. S. (1969). Variance of weighted regression

estimators when sampling errors are independent and heteroscedastic. J. Amer.

Statists. Assoc. 64: 1369-1382.

Cook, D. R. (1977). Detection of influential observations in linear

regression. Technome tries, 19: 15.

Daniel, C; Wood, F. S. (1977). Fitting equations to data. New York: John

Wiley and Sons.

Dennis, J. E. Jr.; Gay, D. M.; Welsch, R. E. (1981a). An adaptive nonlinear

least-squares algorithm. ACM Trans. Math. Software, 7(3): 348-368.

Dennis, J. E. Jr.; Gay, D. M.; Welsch, R. E. (1981b). Algorithm 573:

NL2S0L - an adaptive nonlinear least squares algorithm [E4] . ACM Trans. Math.

Software, 7(3): 369-383.

Dennis, J. E. Jr.; Schnabel, R. B. (1983). Numerical methods for

unconstrained optimization and nonlinear equations. New Jersey:

Prentice-Hall, Inc.

Donaldson, J. R.; Tryon, P. V. (1983). Introduction to STARPAC, the

Standards time series and regression package. Nat. Bur. Stand. (U.S.) Tech.

Note 1068-1.

Dongarra, J. J.; Moler, C. B.; Bunch, J. R.; Stewart, G. W. (1979). UNPACK

users' guide. Philadelphia: SIAM.

Draper, N. R.; Smith, H. (1981). Applied regression analysis, second

edition. New York: John Wiley and Sons.

Filliben, J. J. (1977). User's guide to the DATAPAC data analysis package.

Unpublished.

Fox, P. A.; Hall, A. D.; Schryer, N. L. (1978a). Algorithm 528: framework

for a portable library [Z] . ACM Trans. Math. Software, 4(2): 177-188.

Fox, P. A.; Hall, A. D.; Schryer, N. L. (1978b). The PORT mathematical

subroutine library. ACM Trans. Math. Software, 4(2): 104-126.

Gill, P. E.; Murray, W.; Wright, M. H. (1981). Practical optimization. New

York: Academic Press.

2-R.l
Hoaglin, D. C; Welsch, R. E. (1978). The hat matrix in regression and

ANOVA. American Statistician, 32: 17.

Hogben, D. ; Peavy , S. T.; Varner, R. N. (1971). OMNITAB II user's reference

manual. Nat. Bur. Stand. (U.S.) Tech. Note 552.

Kennedy, W. J. Jr.; Gentle, J. E. (1980). Statistical computing. New York:

Marcel Dekker, Inc.

Lawson, C; Hanson, R.; Kincaid, D. ; Krogh, F. (1979). Basic linear algebra

subprograms for FORTRAN usage. ACM Trans. Math. Software, 5(3): 308-371.

Ryder, G. B. (1974). The PFORT verifier. Software Practice and Experience,

4(4); October-December: 359-377.

Schnabel, R. B. (1982). Finite difference derivatives - theory and practice.

(Unpublished - available from J. R. Donaldson.)

Tryon, P. V.; Donaldson, J. R. (1978). STATLIB: A library of Fortran

subroutines for statistical analysis of experimental data.

(Unpublished - available from J. R. Donaldson.)

2-R.2
NBS-114A rev. 2-aci U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions;

1. PUBLICATION OR

REPORT NO. NBS TN-1068-^

2. Performing Organ. Report No 3. Publ ication Date October 1983

4. TITLE AND SUBTITLE

Nonlinear Least Squares Regression Using STARPAC, The Standards Time Series and

Regression Package

5. AUTHOR(S) Janet R. Donaldson and Peter V. Tryon

6. PERFORMING OT. A n I Z ATI ON (//"joint or other than NBS. see instructions;

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

Report #2 in a series

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City, State, ZIP;

10. SUPPLEMENTARY NOTES

Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a si gnificant

bibliography or literature survey, mention it here)

STARPAC, the Standards Time Series and Regression Package, is a library of Fortran

subroutines for statistical data analysis developed by the Statistical Engineering

Division (SED) of the National Bureau of Standards (NBS), Boulder, Colorado. Earlier

versions of this library were distributed by the SED under the name STATLIB fTryon

and Donaldson, 1978]. STARPAC incorporates many changes to STATLIB, including

additional statistical techniques, improved algorithms and enhanced portability.

STARPAC emphasizes the statistical interpretation of results, and, for this reason,

comprehensive printed reports of auxiliary statistical information, often in

graphical form, are automatically provided to augment the basic statistical

computations performed by each user-callable STARPAC subroutine. STARPAC thus

provides the best features of many stand-alone statistical software programs within

the flexible environment of a subroutine library. This Note documents 16 subroutines

for nonlinear least squares regression. Twelve of these compute the least squares

estimates, performing either weighted or unweighted analysis with either numerically

approximated or user-supplied (analytic) derivatives. The other four are user-callable subroutines for two procedures used within the estimation code: the

first selects optimum step sizes for approximating the partial derviatives of the

model; and the second checks the validity of a user-supplied derivative subroutine.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

derivative checking; NL2SOL; nonlinear least squares; nonlinear regression:

quasi-Newton methods; STARPAC; statistical computing; statistical subroutine library

statistics; derivative step size selection; weighted nonlinear least squares.

13. AVAILABILITY

£g Unlimited

7J For Official Distribution. Do Not Release to NTIS

[X] Order From Superintendent of Documents, U.S^Government Printing Office, Washington, D.C. 20402. n Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

64

15. Price $4.50

fru.S. GOVERNMENT PRINTING OFFICE: 1983-776-876/62 USCOMM-DC 6043-P80



NBS TECHNICAL PUBLICATIONS

PERIODICALS

JOURNAL OF RESEARCH—The Journal of Research of the National Bureau of Standards reports NBS research and develop- ment in those disciplines of the physical and engineering sciences in which the Bureau is active. These include physics, chemistry,

engineering, mathematics, and computer sciences. Papers cover a broad range of subjects, with major emphasis on measurement

methodology and the basic technology underlying standardization. Also included from time to lime are survey articles on topics

closely related to the Bureau's technical and scientific programs. As a special service to subscribers each issue contains complete

citations to all recent Bureau publications in both NBS and non- NBS media. Issued six times a year. Annual subscription: domestic

$18; foreign $22.50. Single copy, $5.50 domestic; $6.90 foreign.

NONPERIODICALS

Monographs—Major contributions to the technical literature on

various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safety codes) developed in cooperation with in- terested industries, professional organizations, and regulatory

bodies. Special Publications—Include proceedings of conferences sponsored by NBS, NBS annual reports, and other special publications

appropriate to this grouping such as wall charts, pocket cards, and

bibliographies

Applied Mathematics Scries—Mathematical tables, manuals, and

studies of special interest to physicists, engineers, chemists,

biologists, mathematicians, computer programmers, and others engaged in scientific and technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical properties of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide program coordinated by NBS under

the .authority of the National Standard Data Act (Public Law

90-396).

NOTE: The principal publication outlet for the foregoing data is the Journal of Physical and Chemical Reference Data (JPCRD)

published quarterly for NBS by the American Chemical Society

(ACS) and the American Institute of Physics (AIP). Subscriptions,

reprints, and supplements available from ACS, I 155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials, components,

systems, and whole structures. The series presents research results,

test methods, and performance criteria related to the structural and

environmental functions and the durability and safety characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of a subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures

published by the Department of Commerce in Part 10, Title 15, of

the Code of Federal Regulations. The standards establish nationally recognized requirements for products, and provide all concerned interests with a basis for common understanding of the

characteristics of the products. NBS administers this program as a supplement to the activities of the private sector standardizing

organizations. Consumer Information Series—Practical information, based on N BS research and experience, covering areas of interest to the con- sumer. Easily understandable language and illustrations provide

useful background knowledge for shopping in today's tech- nological marketplace. Order the above N BS publications from: Superintendent of Documents, Government Priming Office. Washington. DC 20402

Order the following NBS publications— Fl PS and NBSIR's—from

the National Technical Information Service . Springfield. VA 22161

Federal Information Processing Standards Publications (FIPS

PUB)— Publications in this series collectively constitute the Federal Information Processing Standards Register. The Register

serves as the official source of information in the Federal Govern- ment regarding standards issued by NBS pursuant to the Federal Property and Administrative Services Act of 1949 as amended.

Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717(38 FR 12315, dated May II, 1973) and Part 6

of Title 15 CFR (Code of Federal Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or

final reports on work performed by NBS for outside sponsors

(both government and non-government). In general, initial dis- tribution is handled by the sponsor; public distribution is by the National Technical Information Service , Springfield, VA 22161,

in paper copy or microfiche form.
U.S. Department of Commerce

National Bureau of Standards

Washington, DC. 20234

Official Business

Penalty for Private Use $300

POSTAGE AND FEES PAID

U S DEPARTMENT OF COMMERC]

COM-215

SPECIAL FOURTH-CLASS RATI BOOK