A 11 ID 5 b7TD5S

NATL INST OF STANDARDS & TECH R.I.C. 1

NBS

PUBLICATIC

A1 11 02679025

Wolfe, Claire M/An Interactive nonlinear QC100 .U5753 NO.1238 1987_V198 C.1_NBS-P ig OF COMMERCE

tandards

NBS Technical Note 1238

An Interactive Nonlinear Least Squares Program

Claire M. Wolfe, Bert W. Rust, Jeffrey H. Dunn, and Ivor E. Brown

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N$

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS Nl

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS "

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS Nl

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS I IS NBS NBS NBS NBS NBS NBS NBS NBS NBS Nl

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS>

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N

NBS NBS NBS National Bureau ofStandards NBS NBS\

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS Nl

MRS NBS NBS NBS NBS NBS NBS NBS NBS NBS<

?5 NBS NBS NBS NBS NBS NBS NBS NBS Nl

NBS NBS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS Nl

QC

100

I. U 5753

No. 1238

1987

C2
Tm he National Bureau of Standards 1 was established by an act of Congress on March 3, 1901. The Bureau's overall J/^ 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 to assure international competitiveness and leadership of U.S.

industry, science arid technology. NBS work involves development and transfer of measurements, standards and related science and technology, in support of continually improving U.S. productivity, product quality and reliability, innovation

and underlying science and engineering. The Bureau's technical work is performed by the National Measurement

Laboratory, the National Engineering Laboratory, the Institute for Computer Sciences and Technology, and the Institute for Materials Science and Engineering.

The National Measurement Laboratory

Provides the national system of physical and chemical 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; provides advisory and research services to other Government agencies; conducts physical and chemical research; develops,

produces, and distributes Standard Reference Materials; provides

calibration services; and manages the National Standard Reference Data

System. The Laboratory consists of the following centers:

• Basic Standards2 • Radiation Research

• Chemical Physics

• Analytical Chemistry

The National Engineering Laboratory

Provides technology and technical services 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 Engineering3

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 divisions:

Information Systems Engineering

Systems and Software

Technology

Computer Security

Systems and Network

Architecture

Advanced Computer Systems

The Institute for Materials Science and Engineering

Conducts research and provides measurements, data, standards, reference

materials, quantitative understanding and other technical information

fundamental to the processing, structure, properties and performance of

materials; addresses the scientific basis for new advanced materials

technologies; plans research around cross-cutting scientific themes such as nondestructive evaluation and phase diagram development; oversees Bureau-wide technical programs in nuclear reactor radiation research and

nondestructive evaluation; and broadly disseminates generic technical information resulting from its programs. The Institute consists of the following Divisions:

• Ceramics

• Fracture and Deformation3 • Polymers

• Metallurgy

• Reactor Radiation

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address

Gaithersburg, MD 20899.

2Some divisions within the center are located at Boulder, CO 80303.

'Located at Boulder, CO, with some elements at Gaithersburg, MD
Research Information Center

National Bureau of Standards

Gaithersburg, Maryland 20899

QLloo

Mo. t*3$

NBS Technical Note 1238

An Interactive Nonlinear Least Squares Program

Claire M. Wolfe, Bert W. Rust, Jeffrey H. Dunn, and Ivor E. Brown

Center for Applied Mathematics

National Engineering Laboratory

National Bureau of Standards

Gaithersburg, MD 20899

July 1987

V

•"»t*u of *

U.S. Department of Commerce

Malcolm Baldrige, Secretary

National Bureau of Standards

Ernest Ambler, Director
National Bureau of Standards U.S. Government Printing Office For sale by the Superintendent

Technical Note 1238 Washington: 1987 of Documents,

Natl. Bur. Stand. (U.S.), U.S. Government Printing Office,

Tech. Note 1238 Washington, DC 20402

160 pages (July 1987)

CODEN: NBTNAE
Abstract

INVAR is an interactive computer code which uses the Stanford variableseparable nonlinear regression program VARPRO to solve nonlinear (and

linear) least squares problems. The variable-separable feature of VARPRO

makes it attractive to users with real-world fitting problems because it it- erates only on the parameters which appear nonlinearly in the model. Not

only does this simplify the iteratation, but it also means that the user is not required to supply initial estimates for the parameters which appear

linearly.

INVAR implements VARPRO within an enviroment providing the user with on-line feedback and the opportunity to make changes, transformations,

and corrections in real-time. It provides extensive statistical diagnostics

and plots of the results. It can provide publication-quality DISSPLA plots

either at a graphics terminal during exectution or at a hard copy device

afterwards. The code is written mostly in standard Fortran with only a few instructions specific to the NBS Cyber 855. With minor changes it can

easily be implemented on other computers.

This report is both a tutorial guide for beginners and a reference manual for experienced users who wish to make changes in the code. It contains

three completly solved example problems. Three appendices contain information necessary for making changes in the programs.

11

Contents

1 INTRODUCTION 3

2 SOME LEAST SQUARES NOTATION 7

2.1 An Introduction 7

2.2 Writing the Model for INVAR 7

2.3 Weighted Least Squares 9

3 THE USER SUPPLIED SUBROUTINES 10

3.1 General Information 10

3.2 Subroutine CALC 11

3.3 Subroutine DERIV 12

3.4 Subroutine WT 13

3.5 Subroutine TRANS 13

3.6 Subroutine BAKTRN 14

4 THE LOGIC OF THE INVAR DRIVER PROGRAM 15

5 EXAMPLE 1 - LIGHT CURVE OF A SUPERNOVA 18

5.1 The Model and the Data 18

5.2 The User-Supplied Subroutines 19

5.3 Input-Output Options 22

5.4 The Terminal Session 23

6 RUNNING INVAR AT THE NATIONAL BUREAU OF

STANDARDS 39

6.1 The INVOKE,INVAR Procedure on the Cyber 855 39

6.2 Other Ways to Run INVAR 40

1
7 PROMPTS, OPTIONS AND DIAGNOSTIC MESSAGES 42

7.1 General Information 42

7.2 Data Initialization Prompts 42

7.3 The Command Point 46

7.4 Running The Least Squares Program 48

7.5 Running the Statistics Programs 49

7.6 The VARPRO Error Diagnostic Messages 49

8 STATISTICAL CONSIDERATIONS 52

9 A TIME SERIES EXAMPLE - THE WILDLIFE CYCLE 55

9.1 The Model and the Data 55

9.2 The User-Supplied Subroutines 60

9.3 The Terminal Session 63

9.4 The Graphical Output 82

10 AN EXAMPLE WITH TWO INDEPENDENT VARIABLES

- VISCOSITY OF ORGANIC COMPOUNDS 89

10.1 The Model and the Data . 89

10.2 The User-Supplied Subroutines 90

10.3 The Terminal Session 95

10.4 The Graphical Output 115

A THE SUBROUTINES IN THE INVAR CODE 124

B CHANGING THE INVAR PROGRAMS 127

B.l Changing the PROGRAM Statement 127

B.2 Changing the DISSPLA Plotting Device Nomination 127

B.3 Changing the Array Dimensions 128

B.4 Changing the Iteration Specification Parameters 137

C DESCRIPTION OF VARIABLES USED IN INVAR 139
Chapter 1

INTRODUCTION

INVAR is an interactive program which solves linear or nonlinear least

squares problems. Its main advantage over other least squares programs

is that it separates the linear and nonlinear parameters, iterating only on

the nonlinear ones. This eliminates the need for initial guesses of the linear

parameters, and usually improves the convergence properties for nonlinear

problems. INVAR was built around the VARPRO subroutine package written originally by John Bolstad [l] at Stanford University, using an algorithm

developed by Golub and Pereya [7]. So far as we know, VARPRO is the only

widely available, variable-separable nonlinear least squares code. Numerical

analysts and statiticians have been slow to see the great practical advantages of a variable-separable algorithm, but scientists in other fields, with

real world modelling problems, very quickly see the advantage and usually

become enthusiatic VARPRO users after their first encounter with the code.

The original VARPRO was not widely used because it was difficult to implement. The user needed a detailed knowledge of FORTRAN in or- der to write a driver, for the program and to alter the subroutine ADA

which performed calculations specific to the individual problem. In addition, VARPRO was designed to be run in batch mode, involving much time

and effort to change initial estimates of the nonlinear parameters and rerun

the program.

J. E. Rail and R. E. Funderlic, who first implemented INVAR [12],

greatly enhanced the utility of the original VARPRO routines. They wrote

a driver for VARPRO which changed it from a batch to an interactive program. The driver provided greater flexibility to the program, allowing the

user to try many different initial estimates all within one terminal session.
Instead of altering the subroutine ADA to fit the particular problem, Rail

and Funderlic rewrote ADA to call for three very simple subroutines. These

subroutines are written by the user and provide all the necessary problemspecific calculations. These subroutines simplified the coding required by

the user and allowed users with a minimal knowledge of FORTRAN to use the program with ease. Other improvements found in the original INVAR include:

• an option allowing the use of numerical derivatives instead of calculating derivatives in closed form,

• extensive statistical diagnostics which can be applied to the results of

the least squares fit,

• rudimentary graphics.

With excellent documentation on how to use the program [12], Rail and Funderlic developed VARPRO into a useful tool for solving linear and nonlinear

least squares problems.

The present implementation of INVAR has been adapted from Rail and

Funderlic's version. The least squares method used by Bolstad and the

statistics added by Rail and Funderlic have been left virtually intact. One

programming error found in the subroutine STAT has been corrected so that the variable (RMS) (residual mean square) is always assigned a value

for later calculations regardless of which statistical tests the user chooses to

run.This version of INVAR has been developed to handle data more effi- ciently than the original. The program now allows the user many more options for manipulating the data. Among the more important of these options

is the ability to transform both the independent and dependent variables before analysis. The user chooses the desired transformation by writing two

simple subroutines which transform and back-transform the data.

Another major improvement is in the graphical output available to the

user. The program can now be run on graphics terminals, or graphical out- put can be directed to other special plotting devices or to a conventional line printer. Graphs to be plotted on plotting terminals and hardcopy devices

are created using the D1SSPLA graphics package. The plots for CRTs and

line printers are conventional line printer graphs, printed in a horizontal for- mat or in a vertical format for time series data. The subroutines to create these plots have been adapted from those of the STATLIB library of P. V.

Tryon and J. R. Donaldson [13].
This document is intended both as a tutorial guide for users not familiar

with INVAR and as a reference manual for experienced users. New users should read the first 5 chapters before attempting to use the code.

Chapters 2 and 3 give the necessary information to write the user supplied subroutines. Chapter 2 deals with writing the model in variable separable form, while Chapter 3 describes the required FORTRAN subroutines

in detail.

Chapter 4 gives a brief description of the logic of the INVAR driver

program. It should be read by beginning users before trying to execute the

program.

Chapters 5, 9 and 10 give examples of real world modelling problems

solved with INVAR, including output from actual terminal sessions. The

problem in Chapter 5 is simpler than the other two. Chapter 9 describes a time-series modelling problem and Chapter 10 describes a problem with

two independent variables.

Chapter 6 describes some Cyber 855 Procedures that can be used to execute INVAR at the National Bureau of Standards and will probably be

of interest only to NBS users. Chapter 7 gives detailed explanations of the prompts and messages that

the user will encounter during the execution of INVAR. Most users should

be able to get started with INVAR after reading Chapter 5, but might want

to refer to Chapter 7 for more details on specific points. Many users will

probably also refer quite often to Table 7.1 which summarizes all of the

INVAR command codes.

The statistical diagnostics generated by INVAR are described in Chapter

8, which is included for users who are not already familiar with the standard

least squares diagnostic plots.

Users who want to change the INVAR programs to meet their own special needs should refer to Appendices A, B and C. Appendix A gives the

subroutine hierarchy for the whole package, Appendix B tells how to make

several changes that are particularly likely to be desired, and Appendix C

defines the variables used in the various subroutines in the package.

The authors would like to acknowledge their indebtedness to the following people whose previous work formed the foundation for this version of

INVAR:

• G. H. Golub and V. Pereyra, for their work on the variable separable

nonlinear least squares algorithm;

• John Bolstad, for his work on the original batch version of VARPRO;
• J. E. Rail and R. E. Funderlic, for their implementation of the original

interactive INVAR program;

• P. V. Tryon and J. R. Donaldson, for the STATLIB line printer plot

routines.

Thanks also go to Ron Boisvert for his help and suggestions in writing the

NOS Procedures for inplementing INVAR on the Cyber 855, and to David

Tate for numerous suggestions to improve the program. Finally, we would

like to thank Drs. Robert E. Funderlic, Takashi Kashwagi, and Larry Reed

for their careful reviews of this manuscript.
Chapter 2

SOME LEAST SQUARES

NOTATION

This chapter is reprinted with slight alterations from [12, pages 5-9].

2.1 An Introduction

It is assumed that the reader is familiar with the method of least squares. The

purpose here is to familiarize the reader with the notation that will help with

the use of INVAR. For motivation to understand this notation, some readers may

immediately want to try to duplicate the simple example in Chapter 4.

2.2 Writing the Model for INVAR

The method of least squares seeks to minimize the sum of the squares of the differ- ence between the observed value and the value predicted by the model equation at each observation point. (The difference is called the residual.) That is, one seeks

to minimize the residual sum of squares m

ro = £(*-*)' (2-1)

t=i

where the j/j are the observed values and the tjj are the predicted values for observations. Each rji is a function of a vector x of unknown parameters and is evaluated

at the observed value ti of the independent variable t. Least squares problems fall into the categories of linear or nonlinear depending

on how the parameters enter into the model. A model such as

r)(x;t) = ii +x2 t + xs t s + x4 e' (2.2)
would be considered linear since the unknown parameters %i can be written as

coefficients of functions which depend only on the independent variable t. Geometrically, the model is linear in the sense that all first order partial derivatives

with respect to the unknown parameters are independent of the parameter values,

and consequently all higher order derivatives are zero. A model which cannot be

written in this form would be considered nonlinear, e.g.,

j?(z; t) = xi + x2 e~

Xit + i4 sin [xs t) + e~

Xtt

. (2.3)

In equation (2.3), 11, z2 and 14 appear linearly whereas 13, 15 and x$ appear

nonlinearly.

Some recent techniques [7] for solving least squares problems take better advantage of the way the unknown parameters fit into the model. This is done not only by taking advantage of the sum of squares form of the problem but also by

considering the models in separable form, i.e., separating the unknown parameters

into ones which are linear and ones which are nonlinear. Thus for INVAR, we must

write our models in the form

ri(c,a;t) = c^or;*) + ••• + cn$n (a;t) + $n+1 (a;t) , (2.4)

where the linear parameters Cj are coefficients multiplying the functions $j of the

nonlinear parameters a and the independent variables t. The separable form can

be used for considering linear problems (no unknown a) and strictly nonlinear

problems (no unknown c) as well as a combination of the two. To emphasize the

difference between the unknown linear and nonlinear parameters, example (2.3) would be rewritten as

rj(e,a;t) = c x + c 2 e" ' 11 + c3 sin {a2 t) + e~ ait

= c 1 $1 (tt;t) + c 2 $2 (a;<) + c s^s(a;0 + <M«;*) ,

where n = 3,

$1 = 1,1 $2 = e~ ait

> $S = sin (a2 t) , *4 = e~

crjt

For this example: To the INVAR question of how many linear parameters are

there, the user would respond 3; to how many nonlinear, 3; to which terms depend

on nonlinear parameter 1, 2; i.e., $2 depends on t*i. Similarly term 3 depends on

nonlinear parameter 2, and 4 on 3. Of course, the user may want to crudely assume n = in (2.4) so that

t] = $i(a;t) , where

$! = «! + a2 e _Qst + a4 sin (as t) + e' a>t

,
and treat all the parameters nonlinearly. In this case INVAR will not be as efficient and the user will have to supply estimates for all the parameters. When the linear and nonlinear parameters are separated, only the latter require estimates.

2.3 Weighted Least Squares

It is often useful to modify the least squares problem by taking into consideration

the "goodness" of the observations. Thus instead of (2.1) we will consider

m

* = £*?(*-*)'• (2-5)

where Wi is a positive weight assigned to the t'th observation (e.g., Wi = 1/<Tj, where

<?,• is the standard deviation of the uncertainty in t/j). According to the notation

here the FORTRAN array WEIGHT{I) is Wf. The user should note that in the

output table that gives the weighted residual, the quantity

m(yi - rii) = (WEIGHT(I) * *l/2) * (Y(I) - ETA(I))

is printed (cf. the sample problems in Chapters 4 and 7).

It has been the authors' experience that there is considerable confusion as to what a weight is. It is our hope that distinguishing between W and WEIGHT

will decrease this confusion. (We changed the name of W to WEIGHT in several places in the coding to accomplish this consistency.)

Linear least squares problems are often written in the following form: Find x such that

r = mjn {\\Ax-y ||> , (2.6)

where A is the m x n matrix of the functions of the independent variables (i.e., the

coefficients of the unknown parameters) at each point, x is the n-vector of unknown

parameters, y is the m-vector of observed values, and || • • •

|| is the Euclidean norm. Thus to put (2.5) in the notation of (2.6) the problem is to minimize | W[Ax — b) || over all x in R. n

, where

W = diag(W1 ,W2 ,...,Wrn ).
Chapter 3

THE USER SUPPLIED

SUBROUTINES

3.1 General Information

In order to use INVAR the user must supply at least one very simple Fortran

subroutine, called CALC, which computes the functions $ and 77: CALC a subroutine which, for a given data point (U,yi) and the current

values of the nonlinear parameters a, computes the functions 4> and

n- Four additional optional subroutines can also be supplied. These four subroutines are: DERIV a subroutine for calculating the partial derivatives of the $ with

respect to the nonlinear paramenters a, WT a subroutine for calculating or otherwise specifying the weights W{

used to define the weighted sum of squared residuals [cf., Equation

(2-5)],

TRANS a subroutine for making mathematical transformations on the

data before making the fit, BAKTRN a subroutine used in conjuction with TRANS to back-transform

the fit obtained.

These four can be dummy routines if the user so chooses, but it is necessary

to create one file containing all five subroutines. After the subroutines are

10
written, they must be compiled, and loaded with the rest of the INVAR

program.

All of the subroutines except for CALC are optional. However, any

subroutine not used in the problem must be present in the user's element

as a dummy subroutine. A dummy subroutine consists of the statement

SUBROUTINE (along with the required name and argument list), followed

by the statements RETURN and END. For examples of the subroutines see Chapters 5, 9 and 10.

3.2 Subroutine CALC

The subroutine CALC calculates the predicted value of the dependent variable for each observation. To write CALC the user must first write the

model of the problem in separable form (cf. Chapter 2). The vectors ALF

and C contain the values for the nonlinear parameters and linear parameters,

respectively. The array T contains the values of the independent variables,

where T(K,J) is the Jth independent variable at the Kth observation. The

variable K is passed into CALC as a parameter, while J must be specified

in the subroutine. The arguments NL, N, MMAX and IV are all used to dimension the arrays in the subroutine. NL is the number of nonlinear

parameters, N is the number of linear parameters, MMAX is the leading

dimension of the T array and IV is the number of independent variables.

The first four statements in CALC must be written as follows:

SUBROUTINE CALC (K ,NL .N.MMAX , IV, C .ALF ,T ,ETA)

PARAMETER (MAXN1=11)

DIMENSION C(N), ALF(NL) , T(MMAX.IV)

COMMON /BLKPH/ PH(MAXNl)

Next, the separate functions $ of the model are assigned to the vector

PH. All PH(I) which are independent of the nonlinear parameters ALF must

appear first, followed by the remaining PH(I). If the moded contains a function $n+ i which does not have a linear coefficient, then the corresponding

value PHI(N+1) should be be final element in the array. Following the cal11
culation of the separate functions, their linear combination, ETA, should be

assigned a value. The general format of this statement is

ETA = C(1)*PH(1) + C(2)*PH(2) + ... + C(N)*PH(N) + PH(N+1)

See Chapter 2 for instructions on how to construct this equation.

3.3 Subroutine DERIV

The subroutine DERIV need be written only when some or all of the partial derivatives of the PH functions are to be calculated in closed form.

If all derivatives are to be calculated numerically, a dummy subroutine is sufficient. DERIV assigns the partial derivatives to the array DPH, where

DPH(I,J) is the partial derivative of the Ith function PH(I) with respect to the Jth nonlinear parameter ALF(J).

The first eight statements must be written as follows:

SUBROUTINE DERIV (K.T.MMAX , IV.ALF ,NL ,NP1)

PARAMETER (MAXN1=11, MAXKG=8)

DIMENSION T(MMAX, IV) , ALF(NL)

COMMON /BLKDPH/ DPH(MAXN1 .MAXKG)

DO 10 J=1,NL

DO 10 I=1,NP1

DPH(I.J) = 0.0

10 CONTINUE .

The DO loop serves to initialize all the partial derivatives to zero, and is not

necessary when using a dummy routine. Parameters K, MMAX, IV, and

NL are the same as in subroutine CALC, and NP1 = N + 1. The remaining

12
statements depend upon the specific problem, but each derivative the user wishes to calculate in closed form must be assigned a value here.

3.4 Subroutine WT

The subroutine WT, also optional, calculates the weights for the problem.

If the user does not wish to calculate weights, a dummy subroutine must

still be provided. In this case, weights can still be entered with the data

points, or an equal weighting ( WEIGHT(I) = 1) will be assumed. In WT,

the vector Y contains all the observed values of the dependent variable. WT

should assign the weighting values to the vector WEIGHT.

The first two statements of WT must be written as follows:

SUBROUTINE WT (M, WEIGHT, Y)

DIMENSION WEIGHT(M), Y(M)

The remaining statements may be written by the user to suit the individual

problem, but a weight must be assigned for each of the M observations of

Y.

3.5 Subroutine TRANS

If the user wants to transform the independent and/or the dependent variables before the fit, the subroutine TRANS must be written. Otherwise, a dummy subroutine will suffice. The T array and the Y vector contain the in- dependent and the dependent variables, respectively. On entering TRANS,

these arrays contain the original data, and on exiting TRANS, they contain

the transformed data.

The first two statements must be written as follows:

SUBROUTINE TRANS (T.BT.Y ,BY , IV.M.MMAX)

DIMENSION T(MMAX.IV), BT(MMAX.IV), Y(M). BY(M)

13
TRANS should not only transform variables but should also retain their

original values (the "back-transformed" values) in array BT for the independent variable and vector BY for the dependent variable. Therefore, the

user must include DO loops containing the statements

BT(I.J) = T(I,J)

BY(I) = Y(I)

before the values are transformed. Even if only one of the variables is to be

transformed, the user must always assign these "back-transformed" values to both variables. The user then transforms the variables as desired, assigning

the transformed values to the original T array and Y vector.

3.6 Subroutine BAKTRN

The subroutine BAKTRN must also be written if transformed values are

desired. Otherwise, BAKTRN need only be a dummy subroutine. If data

has been transformed, the predicted values in vector ETA are calculated in

CALC using the transformed data. INVAR plots the observed and predicted

values twice, once as transformed data, and once as back-transformed data.

Therefore, ETA must be back-transformed in order to be plotted with the

original values of the dependent variable. BAKTRN performs this backtransformation of ETA, placing the new values in the vector BETA.

The first two statements must be written as follows:

SUBROUTINE BAKTRN (ETA, BETA, M)

DIMENSION ETA(M), BETA(M) .

The remainder of the subroutine is written by the user to do the required

back-transformation of ETA. ETA must be backtransformed by exactly the

inverse method used to transform the vector Y in the subroutine TRANS.

14
Chapter 4

THE LOGIC OF THE

INVAR DRIVER

PROGRAM

No attempt is made here to describe the program package in detail. Anyone

interested in these details is referred to Rail and Funderlic's documentation

on INVAR [12]. However, a general understanding of how the program works

is desirable to ensure efficient use of INVAR.

The program's driver is divided into two sections. The first of these

sections is the Initialization Segment which asks for the data and initial values needed to run the least squares program. These values include the

number of data points, linear and nonlinear parameters and independent

variables, along with the data itself (independent and dependent variables

and weights) and initial estimates for the nonlinear parameters. Also the

program must be told which terms in the model depend on each of the

nonlinear parameters, and it must be given increments for the numerical

derivatives if the latter are used.

The second section of the driver is the Command Segment which contains

the bulk of the program. The user enters into this segment at the Command

Point when the program cues "NEXT STEP?". There are four types of

steps which can be taken: (1) print current values of the data or turn on/off

printing flags; (2) change the current parameter estimates, transform the

data, or add/delete data points; (3) run the least squares program or the

statistical analysis; and (4) save the current data in a file and/or exit the

program.

15
All of these steps terminate by returning to the Command Point, except

for those steps which exit the program. It is this looping in the Command

Segment that makes INVAR such a flexible program. The user can run the

least squares program as many times as desired in one execution of INVAR,

each time changing various values and running statistical analyses, in order

to find the best fit. An extremely abbreviated flowchart depicting INVAR's driver is given

in Figure 4.1. Except for the data initialization, all actions of INVAR are

initiated by typing in a command code at the Command Point. Table 7.1 gives a list of all the command codes which may be chosen by the user

at the Command Point. The same table is available interactively during

execution of the program by typing "H" at the Command Point and, after

the program's response, entering "A".

16
(jtST)

Print

options

Chanfi

options

Initialization segment

Figure 4.1: Abbreviated flowchart of the INVAR driver program

17
Chapter 5

EXAMPLE 1 - LIGHT

CURVE OF A

SUPERNOVA

5.1 The Model and the Data

Type I supernovae have been found to follow a specific pattern of luminosity.

The latter part of this pattern may be described by the model

y{t) = Cx exp( ) + C2 exp( ) , Oil <*2

where t is the time (in days) after maximum luminosity and y is the relative

luminosity (relative to max. luminosity). In terms of the notation discussed

in Chapter 2, this model is written as

r/(d, C2 , <*i, «2 ; t) = Ci*i(ai; t) + C2 $2 (a2 ; .

where

$i(ai;t) = exp(-J_) ,

<f> 2 {a2 ;t) = exp{-±).

The model above will be fitted to data collected in 1939 for the type I supernova SN1939A. The data are given in Table 5.1. The peak luminosity

occurred at day 0.0, but all measurements taken before day 7.0 will not be

used here since the above model cannot account for the luminosity before

and immediately after the peak. These data were written, as two columns

18
Table 5.1: Luminosity of SN1939A (Relative to Peak Luminosity)

Day Luminosity Day Luminosity Day Luminosity

7.0 0.6310 19.0 0.1318 57.0 0.05754

7.0 0.8318 20.9 0.1585 85.0 0.03631

14.8 0.2754 25.8 0.1096 109.0 0.02291

16.0 0.1445 26.8 0.1445 110.0 0.02291

16.9 0.2089 28.0 0.09120 141.0 0.01738

17.0 0.1585 53.0 0.06310 142.0 0.01585

18.8 0.1585 54.0 0.06918 168.0 0.009120

of numbers, into a file named SNDATA which will be used in the following

as an input file for the INVAR program.

5.2 The User-Supplied Subroutines

The user-supplied subroutines reqired to make the fit are given in the following paragraph. These subroutines were all written into a single file, called

SNCALC, in accordance with the directions given in Chapter 3. Subroutine

CALC uses the values of PH and ETA obtained from writing the model in separable form. In this example, there is only one independent variable,

t = time (days). DERIV calculates the derivatives from the closed form

expressions

(af?W-i)

|*2 =0 dc*2 = (AW-i).

assigning non-zero values only to DPH(l,l) and DPH(2,2) since PH(l) depends only on ALF(l) and PH(2) depends only on ALF(2). Subroutine WT

assigns a weighting of l/yf to each data point (ti, yt ) in order to assure that

the observations make roughly equal contributions to the determination of

the fitting parameters. This is a very common procedure when the y-values

span several orders of magnitude.

The subroutines for this example are:

19
c

c

c

c

c

c

c

c

c

c-c

c

c

c

c

c

c

SUBROUTINE CALC ( K, NL. N. MMAX . IV. C. ALF, T. ETA )

PARAMETER ( MAXN1 = 11 )

DIMENSION C(N). ALF(NL) . T(MMAX.IV)

COMMON /BLKPH/ PH(MAXNl)

PH(1) = EXP(-T(K.l) / ALF(l))

PH(2) = EXP(-T(K,1) / ALF(2))

ETA = C(1)*PH(1) + C(2)*PH(2)

RETURN

END

SUBROUTINE DERIV ( K. T, MMAX. IV, ALF. NL, NP1 )

PARAMETER ( MAXN1 = 11 , MAXKG = 8 )

DIMENSION ALF(NL).T(MMAX,IV)

COMMON /BLKDPH/ DPH(MAXN1 , MAXKG)

DO 10 J=1,NL

DO 10 1=1, NP1

DPH(I.J) = 0.0

10 CONTINUE

C

C DPH(I.J) SHOULD BE SET TO THE PARTIAL DERIVATIVE

C OF PH(I) WITH RESPECT TO ALF(J)

C

20
DPH(l.l) = T(K,1) / (ALF(1)**2) * EXP(-T(K.1)/ALF(1))

DPH(2.2) = T(K.l) / (ALF(2)**2) * EXP(-T(K,1)/ALF(2)

)

C

RETURN

END

C

C

C

- cc SUBROUTINE WT ( M. WEIGHT, Y ) C

C

DIMENSION Y(M).WEIGHT(M)

C

DO 10 1=1,

N

WEIGHT(I) = 1.0/Y(I)**2

10 CONTINUE

C

RETURN

END

C

C

C ---

c

c SUBROUTINE TRANS ( T, BT, Y, BY. IV, M, MMAX ) C

C

DIMENSION T(MMAX, IV), BT(MMAX.IV), Y(M) , BY(M)

C

C-- DUMMY SUBROUTINE

C

RETURN

END

C

C

C

c

21
c SUBROUTINE BAKTRN ( ETA, BETA, M ) C

C

DIMENSION ETA(M), BETA(M)

Cc DUMMY SUBROUTINE

C

RETURN

END

No transformations of the data are required so the subroutines TRANS and

BAKTRN are dummy subroutines.

The above subroutines were written into a file named SNCALC which was compiled and loaded with the INVAR programs and then executed in a terminal session to be described in the following sections.

5.3 Input-Output Options

There are basically three options for inputting data to INVAR:

1. Type all of the data, including the observations (£,-, y;), into the program during execution,

2. Type in only the specifications for the 3> functions and their derivatives

with respect to the a parameters, but read the observations (U,yi)

from a data file,

3. Read all of the specifications and the observations from a file that was

previously written by the INVAR program.

Initially, one of the first two options must be used, but once this has been

done, INVAR allows the user to save the resulting data set in a special

"VARPRO defined" file so that the final option can be used in subsequent

runs with the same data set. This third option reduces the amount of typing

for the user. The second option was used in the terminal session described

in the next section.

Because INVAR is an interactive program the basic output file is the

user's terminal, but the program also provides the option of writing all the

results to a file (unit 10) which is saved for later processing. If the user

22
exercises this option he can obtain a hard copy of his results on the printer

of his choice.

There are actually two versions of INVAR: INVAR2 which uses the DISSPLA plotting package and INVARl which provides only line-printer graphs

of the results. INVAR2 also provides line printer graphs in case the program crashes and the DISSPLA plots are lost. If INVAR2 is chosen, the

user has two ways to get the graphs. If he is running on a Tektronics graphical terminal he can receive the graphs interactively during execution of the

program. In this case the program uses CALL PTEKAL to initialize the DISSPLA plotting package. The other option is to write the plots to a file which

will be plotted on a QMS laser printer after the terminal session has been

completed. In this case the program uses CALL TK4014 and CALL HWSPEC

to initialize DISSPLA. Users can easily change the program to accomodate

other plotting devices. Instructions for making the required changes are given in Appendix B.

5.4 The Terminal Session

The sample terminal session output that follows was obtained by executing

INVAR2 on the NBS CYBER 855. It illustrates only some of the options

available in INVAR. (For a more comprehensive example, see Chapter 9.)

The program was run from a Tektronics plotting terminal, using the option

of interactive DISSPLA plots and requesting line printer output as well as the normal terminal output. Before reading this example, the following

points should be noted:

1. The first 17 non-blank lines have nothing to do with the INVAR program. They were generated by the NOS procedures used to initialize execution on the NBS Cyber 855. The first line generated by INVAR

is the command "ENTER TITLE (MAXIMUM OF 40 CHARACTERS):".

2. The question mark at the beginning of each line immediately following an INVAR request for information is the Cyber 855 prompt to remind

the user to type in something.

3. Prompts requiring a yes or no response were answered with a "Y" for yes and a "N" for no, although any character other than "Y" will be

interpreted as no.

4. The four queries immediately following the choice of "PLOTS SHOWN

INTERACTIVELY ON A PLOTTING SCREEN" were issued by the Tectron23
ics software called by the DISSPLA initialization "CALL PTEKAL".

5. The 21 observations (U,yi) were read from the file SNDATA.

6. The program requires initial estimates for all of the nonlinear parameters. The initial estimates used for cti and a2 were 5 and 60 days,

respectively.

7. When the model and the data have been specified, INVAR prints the

prompt "NEXT STEP? (TYPE "H" FOR HELP)" and waits for the user

to type in a command code for one of several possible tasks. When

such a task has been completed, INVAR will again print the above

prompt for the next task. A complete list of the possible tasks and

their command codes is given in Table 7.1. The present example used

the command codes "RV, for running the VARPRO least squares

programs, "RS", for running a statistical analysis of the results, and

"SO", to stop without saving the input data in a VARPRO defined file.

8. The Tektronics software clears the screen before making each plot and

then the program stops awaiting a response from the user. The Cyber

855 prints a question mark in the upper left corner of the screen to indicate that some response is expected. In the present example, we

responded with a "C followed by a carriage return because a carriage

return by itself is interpreted by the 855 as an end of file, causing the

program to crash.

9. When using interactive Tektronics plotting, it is a good idea for the

user to clear the screen after each plot since otherwise new INVAR

prompts or results may be written over the preceding plot.

Chapters 6 and 7 give more detailed explanations of the prompts and

responses. The reader may find it useful to refer to them while reading the

following sample terminal session.

/INVOKE, INVAR

INTERACTIVE NONLINEAR LEAST SQUARES PROGRAM

PROC,INVARX,CALCF[,qMSADDR][,DATAF] [.INV102] .

Executes one of two versions, INVAR1 or INVAR2. of INVAR, an interactive, variable separable, nonlinear least squares program.

24
Local file CALCF contains FORTRAN coded, user-supplied subroutines

Hardcopy output and plots can be sent to QMS laser printer with

address QMSADDR.

Permanent file DATAF contains data to be fit (not needed if data

already in a local file, or is to be typed in during exection)

.

Version to be executed specified by parameter INV102. Default is

INVAR2 which doeB DISSPLA plots. INVAR1 does not use DISSPLA.

PARAMETERS FOR INVARX ARE CALCF, QMSADDR, DATAF, INV102

Enter CALCF FORTRAN file of user-supplied subroutines? SNCALC

Enter QMSADDR QMS Address for output (Default = NONE)? B225

Enter DATAF Data file (Default = NONE)? SNDATA

Enter INV102 INVAR1 or INVAR2 (Default = INVAR2)?

ENTER TITLE (MAXIMUM OF 40 CHARACTERS)

? LIGHT CURVE FOR SN1939A

DO YOU WANT DISSPLA PLOTS? (TYPE "Y" FOR YES)

? Y

DO YOU WANT THESE PLOTS SHOWN INTERACTIVELY ON A PLOTTING SCREEN (TYPE "1"); OR SAVED IN A FILE TO BE PRINTED OUT LATER (TYPE "2")?

? 1

ENTER DESIRED MODEL

? 4014

ENTER LINESPEED (CHARACTERS PER SECOND)

? 960

ENTER RESOLUTION MODE (0=LOW)

? 1

ENTER OPTION (0=DEF,31=8 PENS ON 4662 , 36=PAPER FEED ON 4663)

?

ENTER LABEL FOR INDEPENDENT VARIABLE 1 (MAX. OF 47 CHARACTERS):

? DAY

ENTER LABEL FOR DEPENDENT VARIABLE (MAX. OF 47 CHARACTERS):

? RELATIVE LUMINOSITY

25
DO YOU WANT LINE PRINTER OUTPUT?

? Y

ARE YOUR DATA, WEIGHTS, INITIAL VALUES, ETC. ALREADY IN A VARPRO -DEFINED FILE ?

? N

NO. OF DATA POINTS:

? 21

NO. OF LINEAR PARAMETERS:

? 2

NO. OF NONLINEAR PARAMETERS:

? 2

NO. OF INDEPENDENT VARIABLES:

? 1

ARE DATA IN A FILE?

? Y

WHAT IS THE FILENAME?- -NO MORE THAN 26 CHARS

(TYPE "QUIT" TO EXIT)

? SNDATA

ARE WEIGHTS TO BE CALCULATED IN SUBROUTINE WT?

? Y

DO DEPENDENT OR INDEPENDENT VARIABLES

NEED TO BE TRANSFORMED BEFORE ANALYSIS?

? N

INITIAL VALUES OF NONLINEAR PARAMETERS

1 ? 5.0

2

? 60.0

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 1: ? 1

WHICH TERMS ARE THEY?

26
? 1

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 2: ? 1

WHICH TERMS ARE THEY?

? 2

ALL DERIVATIVES IN CLOSED FORM?

? Y

WOULD YOU LIKE TO SAVE THE DATA JUST INPUT?

? N

DOES DATA FOLLOW AN EVENLY SPACED TIME SERIES?

? N

NOTE: THE COMMAND "SC" WILL SAVE YOUR DATA FILE AT ANY TIME

NEXT STEP? (TYPE "H" FOR HELP)

? RV

DO YOU WISH TO CHANGE ANY DEFAULTS IN THE LEAST SQUARES PROGRAM?

(50: MAX. NO. OF ITERATIONS)

(l.OE-5: CONVERG. TOLERANCE)

(1: OUTPUT PRINTING INDEX)

? N

1 LIGHT CURVE FOR SUPERNOVA 1939A

WEIGHTED NORM OF RESIDUAL = 6.7995272E-01

NU = .1000000E+01

ITERATION 1 NONLINEAR PARAMETERS

4.B931719E+00 6.0B32666E+01

1 WEIGHTED NORM OF RESIDUAL = 6.4835577E-01

NU = .6000000E+00

NORM(DELTA-ALF) / NORM(ALF) = 1 . 104E-02

27
ITERATION 2 NONLINEAR PARAMETERS

4.2716641E+00 6.0017367E+01

1 WEIGHTED NORM OF RESIDUAL = 6.3826BB7E-01

NU = .2600000E+00

NORM(DELTA-ALF) / NORM(ALF) = 1.010E-02

ITERATION 3 NONLINEAR PARAMETERS

4.1469268E+00 5 .9472842E+01

1 WEIGHTED NORM OF RESIDUAL = 6.3711B32E-01

NU = .12BOOOOE+00

NORM (DELTA- ALF) / NORM(ALF) = 9.370E-03

ITERATION 4 NONLINEAR PARAMETERS

4.1228B19E+00 6 .93B3197E+01

1 WEIGHTED NORM OF RESIDUAL = 6. 3707619E-01

NU = 6250000E-01

NORM(DELTA-ALF) / NORM(ALF) = 2.0B1E-03

ITERATION B NONLINEAR PARAMETERS

4.1202209E+00 6 .93398B8E+01

1 WEIGHTED NORM OF RESIDUAL = 6.3707674E-01

NU = .312B000E-01

NORM(DELTA-ALF) / NORM(ALF) = 2 . 286E-04

ITERATION 6 NONLINEAR PARAMETERS

4.1199966E+00 B . 9338726E+01

1 WEIGHTED NORM OF RESIDUAL = 6 . 3707673E-01

NU = .1B62B00E-01

NORM(DELTA-ALF) / NORM(ALF) = 1.941E-0B

ITERATION 7 NONLINEAR PARAMETERS

28
4.1199788E+00 6 . 9338638E+01

1 WEIGHTED NORM OF RESIDUAL = 6.3707573E-01

NU = .7812500E-02

NORM(DELTA-ALF) / NORM(ALF) = 1.B17E-06

LINEAR PARAMETERS

.3069174E+01 . 1B86064E+00

NONLINEAR PARAMETERS

.4119979E+01 6933864E+02

WEIGHTED NORM OF RESIDUAL = .6370767E+00

WEIGHTED ESTIMATED VARIANCE = .2387444E-01

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED

BY TYPING "Y": TABLE OF OBSERVED ft PREDICTED VALUES, WTD RESIDUAL:

? Y

COVARIANCE MATRIX:

? Y

CORRELATION MATRIX:

? Y

COVARIANCE MATRIX

ROW COLUMN 1 COLUMN 2 COLUMN 3 COLUMN 4

1 9.1B6764E-01 7.486032E-03 -B .032086E-01 -1 .393200E+00

2 7.486032E-03 2.77B18BE-04 -6. B79036E-03 -6 . B168B6E-02

3 -6.032086E-01 -6 . B79036E-03 3.362799E-01 1.23B1B4E+00

4 -1.393200E+00 -B .B168B6E-02 1.23B1B4E+00 1.371202E+01

1

INDEX IND VAR 1 OBS VAL PRED VAL WTD RESID WEIGHTS

29
1 7 00000E+00 6 31000E-01 7 0218BE-01 -1. 12813E-01 2 B11B5E+00

2 7 00000E+00 8 31800E-01 7 0218BE-01 1. 6B826E-01 1. 44B31E+00

3 1. 48000E+01 2 7B400E-01 2 08109E-01 2 44341E-01 1, 31848E+01

4 1 60000E+01 1 44500E-01 1 84279E-01 -2 75288E-01 4 78921E+01

5 1 69000E+01 2 08900E-01 1 70062E-01 1. 86916E-01 2 29152E+01

6 1 70000E+01 1 68600E-01 1 68644E-01 -6 39993E-02 3 98064E+01

7 1 88000E+01 1 68600E-01 1 47B48E-01 6 91004E-02 3 98064E+01

8 1 90000E+01 1 31800E-01 1 .46642E-01 -1 0B023E-01 B 7B664E+01

9 2 .09000E+01 1 68600E-01 1 30748E-01 1 7B093E-01 3 980B4E+01

10 2 68000E+01 1 09600E-01 1 08B35E-01 9 718B0E-03 8 32490E+01

11 2 68000E+01 1 44600E-01 1 0BBB7E-01 2 69498E-01 4 78921E+01

12 2 80000E+01 9 . 12000E-02 1 .02376E-01 -1 22B43E-01 1 20229E+02

13 5 .30000E+01 6 31000E-02 6 .49339E-02 -2 90630E-02 2 .B115BE+02

14 5 .40000E+01 6 .91800E-02 6 .38472E-02 7 70863E-02 2 08948E+02

15 6 .70000E+01 5 .76400E-02 6 06966E-02 -5 48B84E-02 3 . 02037E+02

16 8 .B0000E+01 3 .63100E-02 3 .78628E-02 -4 . 276BBE-02 7 B8486E+02

17 1 . 09000E+02 2 29100E-02 2 B2673E-02 -1 02894E-01 1 90B24E+03

18 1 . 10000E+02 2 .29100E-02 2 .484B1E-02 -8 .44634E-02 1 .90524E+03

19 1 .41000E+02 1 .73800E-02 1 .473B1E-02 1 B2180E-01 3 310BBE+03

20 1 . 42000E+02 1 .B8B00E-02 1 .44889E-02 8 B87B4E-02 3 980B4E+03

21 1 .68000E+02 9 . 12000E-03 9 348B2E-03 -2 . B0B70E-02 1 . 20229E+04

CORRELATION MATRIX

ROW COLUMN COLUMN COLUMN COLUMN

1 1 OOOOOOE+OO 4 69B444E-01

2 4 69B444E-01 1 . OOOOOOE+OO

3 -9 068314E-01 -6 810293E-01

4 -3 931802E-01 -8 .943241E-01

-9.068314E-01 -3.931802E-01

-6.810293E-01 -8 .943241E-01

1. OOOOOOE+OO B.7B2013E-01

B.7B2013E-01 1. OOOOOOE+OO

PLOT ROUTINE CURRENTLY SET NOT TO DO ANY PLOTS.

DO YOU WANT TO START PLOTTING?

? Y

IF YOU WANT TO PLOT RESIDUALS, TYPE "1"

IF YOU WANT TO PLOT OBSERVED ft PREDICTED VALUES. TYPE "2"

IF YOU WANT TO PLOT BOTH, TYPE "3":

? 3

30
7 C

RESIDURLS

a

a

a

a: «$>

a

a

a fc % •

a

a

*' *' N^' $ s*'

N#'s<f>*

t

OTY

i

i

31
7 C

LIGHT CURVE FOR SN1939R

% - t—

tno

_>UJ

l- *'

<r

_jUJ

DC

*• *' <>' <> s

<*> f s# s

<»- N^'

DAY

32
NEXT STEP? (TYPE "H" FOR HELP)

? RS

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED: PARAMETERS, THEIR STD DEV * T-RATIO:

? Y

RES SUM OF SqS, RES MEAN SQ, RES STD DEV ft COEFF OF DETERMINATION:

? Y

TABLE OF OBS. PRED VAL, RESID t STD RESID:

? Y

PLOT OF STD RES VS IND VAR:

? Y

PLOT OF STD RES VS PRED VAL:

? Y

NORMAL PROBABILITY PLOT:

? Y

1 PARAMETER S.D. OF PARM. T-RATIO

3.069174E+00 9.569098E-01 3.207380E+00

1.586064E-01 1.666889E-02 9.B20828E+00

4.119979E+00 B.798965E-01 7.104680E+00

5.933864E+01 3.7O2975E+00 1.602459E+01

WTD RESIDUAL SUM OF SQUARES: 4.058655E-01

WTD RESIDUAL MEAN SQUARE: 2.387444E-02

WTD RESIDUAL STANDARD ERROR: 1.B45136E-01

COEFFICIENT OF DETERMINATION (R- SQUARE) : 9.673397E-01

33
INDEX IND VAR 1 OBS VAL PRED VAL RES ID STD RES

123

4

56

7

89

10

111213

14

15

17

181920

21

16 8

OOOOOE+00

OOOOOE+OO

48000E+01

60000E+01

69000E+01

70000E+01

88000E+01

90000E+01

09000E+01

58000E+01

68000E+01

80000E+01

30000E+01

40000E+01

70000E+01

.60000E+01

. 09000E+02

. 10000E+02

.41000E+02

.42000E+02

.68000E+02

6.31000E-01

8.31800E-01

2.7B400E-01

1.44B00E-01

2.08900E-01

1.68600E-01

1.68B00E-01

1.31800E-01

1.B8B00E-01

1.09600E-01

1.44B00E-01

9.12000E-02

6.31000E-02

6.91800E-02

6.76400E-02

3.63100E-02

2.29100E-02

2.29100E-02

1.73800E-02

1.B8BOOE-02

9.12000E-03

7.0218BE-01

7.0218BE-01

2.08109E-01

1.84279E-01

1.70062E-01

1.68644E-01

1.47B48E-01

1.4B642E-01

1.30748E-01

1.08B35E-01

1.0B5B7E-01

1.02376E-01

6.49339E-02

6.38472E-02

6.06966E-02

78628E-02

B2673E-02

48451E-02

47351E-02

1.44889E-02

9.34852E-03

7.11847E-02

1.2961BE-01

6.7291BE-02

3.97792E-02

3.88379E-02

1.01439E-02

1.09B24E-02

1.38421E-02

2.77B22E-02

1.06B1BE-03

3.8942BE-02

1.11759E-02

1.83388E-03

B.33283E-03

1B6B6E-03

56282E-03

3B730E-03

93506E-03

2.64488E-03

1.36112E-03

2.28B20E-04

-1.54296E-01

6.69810E-01

1 . 89306E+00

-3.91233E+00

1 . 46276E+00

-6.20178E-01

6.16663E-01

-1.20026E+00

1.52620E+00

1.94168E-01

2.48441E+00

- 1 . 60334E+00

-1.0B740E+00

1.03859E+00

-1.96701E+00

-4.61768E-01

-9.36326E-01

-7.61009E-01

1.08B20E+00

6.22178E-01

-1.80473E-01

ENTER ANY INTEGER TO CONTINUE

? 1

34
7 C

STAND. RESIDURLS VS. IND. VARIRBLE

o

- <n v

_i

<T2 N

' inU(K o -

UJ5 S?

a:

'

<EO

<r v

(n

a

o0°

a a a a a o a

a

a

o

o

V *• *' * N

<*>- s

*>- s

rf>- s

«>- f

OflY

>

35
7 C

STRND. RESIDURLS VS. PRED. VRLUES

->*

'

o

- i

mUJ

or o

UJM

O

«» e b

o a

B

B

D

B °

° O B

B

B

8

f a

'** * N *V ^ * «> %

*• ©>* *• ** ** V ** RELATIVE LUMINOSITY

A %

36
7 C

NQRMRL PROBABILITY PLOT

S e O " in~ V n O 3*

D

=> o - o Q s* o'o o

inu

o o Q o UJ _ M

^ ^ " ° o s- o <£'

<r o o

£<v> I-'*

(T) Q " ,V

>*,<&*> O <o O *» O

,V ,N* ,N" ,fc* ^' O* V N" V

EXPECTED STflNOflROIZED RESIDURLS

1

37
NEXT STEP? (TYPE "H" FOR HELP)

? SO

END DF DISSPLA 9.0 -- 6447 VECTORS IN 5 PLOTS.

RUN ON 4/13/87 USING SERIAL NUMBER AT NBS/CS**2 - GAITHERSBURG. MD. PROPRIETARY SOFTWARE PRODUCT OF ISSCO, SAN DIEGO. CA. 2870 VIRTUAL STORAGE REFERENCES; 4 READS; WRITES.

PROCESSING FILE - ZZLPOUT PLEASE WAIT

EXIT.

/

Note that after the user typed in the command code "SO" the INVAR

program was terminated. The next five line were written by the DISSPLA

programs and the line beginning with "PROCESSING FILE" was written by

the QMS procedure for printing the user's hard copy output.

38
Chapter 6

RUNNING INVAR AT

THE NATIONAL BUREAU

OF STANDARDS

6.1 The INVOKE,INVAR Procedure on the Cyber 855

Users at the National Bureau of Standards can run INVAR on the Cyber

855 simply by typing in the command

INVOKE, INVAR

.

This initiates an interactive NOS procedure called INVARX which asks the

user to provide one required and three optional pieces of information.

The required information is the name of the file containing the user- supplied subroutines CALC, DERIV, WT, TRANS, and BAKTRN (see

Chapter 3). If this is not already a local file, INVARX assumes that it

is a permanent file and attempts to GET it for the user. When loading

the INVAR object code with the user-supplied subroutines, the INVOKE

procedure attaches the CMLIB, IMSL and NAG libraries. Thus the user is free to call any routine in these libraries from his own subroutines.

The three optional pieces of information requested by INVARX are:

1. The address of a QMS laser printer for printing the optional hardcopy

(line-printer) output and the optional DISSPLA plots. If an address is specified, the user will always receive cross-referenced Fortran compilation listings of the user-supplied subroutines. The INVAR program

39
will also ask if he wants hardcopy output, and, if the reply is affirmative, the INVAR output will follow the compilation listings. Some users may wish to delay the decision about whether to get hardcopy output

until after the terminal session has been completed. Others may wish

to receive hardcopy on some printer other than a QMS printer. Such

users should not specify a QMS address to the INVARX procedure,

but should answer the INVAR query in the affirmative. When the

program execution is completed, they will receive the message

YOUR HARDCOPY OUTPUT IS IN FILE ZZLPOUT . They can then either discard the local file ZZLPOUT or send it to the printer of their choice. Note, however, if the user exercises the

DISSPLA plotting option (see below) for QMS plots, he must specify a QMS address.

2. The name of the file containing the observational data. If the data file

is already a local file, there is no need to specify anything here. When

the INVAR program asks for the name of the data file the user simply

types it in. If the data are in a permanent file but not in a local file, the user should type the permanent file name so that the INVARX

procedure can GET it for him. The INVAR program will still ask for

the data file name, and the user must type it in again. Recall that

there are two kinds of data files that might be used (see Chapter 5).

3. The version of INVAR to use . There are two versions of INVAR available to the user. INVARl does not have the option of producing DISSPLA graphics while INVAR2 does have that option. If neither version

is specified, INVAR2 is taken as the default. If INVAR2 is used, there

are two ways to receive the DISSPLA plots: on a Tektronics graphical

terminal, or as QMS laser printer plots. If the latter option is chosen

a QMS printer address must be specified as described above.

For an example of the use of the INVOKE,INVAR procedure while in line mode (eg., on a Tektronics terminal) see the sample terminal session

given in Chapter 5. For a example of its use in screen mode see the sample

terminal session in Chapter 9.

6.2 Other Ways to Run INVAR

Many NBS users may not be satisfied with the options provided by the

INVOKE procedure or by the INVAR program, and others may not want

40
to run on the Cyber 855.

The INVAR source code can be obtained by any 855 user by typing

GET, SOURCE/UN = INVAR .

The file SOURCE contains two records: INVARl and INVAR2. Instruc- tions for altering these codes are given in Appendix B. Possible alterations

include:

1. Changing the DISSPLA plotting device nominations,

2. Changing the maximum number of observations,

3. Changing the maximum number of independent variables,

4. Changing the maximum number of linear and/or nonlinear parameters,

5. Changing the iteration specification parameters.

Users who are satisfied with the INVAR codes and are content to run on the 855, but do not like the INVOKE procedure, can write their own

procedures to use one of the object code files

0BJECT1/UN = INVAR or 0BJECT2/UN = INVAR ,

which contain INVARl and INVAR2.

Users who want develop their own procedures by altering the ones used

by INVOKE,INVAR can get the latter from

INVARX/UN = INVAR and IOPROCS/UN = INVAR .

41
Chapter 7

PROMPTS, OPTIONS

AND DIAGNOSTIC

MESSAGES

7.1 General Information

The following list, which will continue through several sections, contains

prompts that INVAR may give during a terminal session. Following each

prompt (which is written in upper case letters) is an explanation of the

option and the various ways in which the user can respond. The prompts

are listed in roughly the same order in which they will be encountered, but

INVAR follows many possible paths of logic and the route taken depends on how the user answers these prompts. Not all prompts will always be given

in one terminal session.

7.2 Data Initialization Prompts

1) ENTER TITLE: (MAXIMUM OF 40 CHARACTERS) The user enters a title for the problem (up to 40 characters long). This title becomes the heading for each run of the least squares program and

the title for the DISSPLA graph of observed and predicted values (if the DISSPLA plotting option is chosen). If INVAR1 is being run, the

program advances to prompt 6). INVAR2 advances to prompt 2).

(Chapter 6 explains the difference between INVAR1 and INVAR2.)

42
2) DO YOU WANT DISSPLA PLOTS? (TYPE KY" FOR YES)

If the answer is not "Y", the program advances to prompt 6).

3) DO YOU WANT THESE PLOTS SHOWN INTERACTIVELY

ON A PLOTTING SCREEN (TYPE "l"); OR SAVED IN

A FILE TO BE PRINTED OUT LATER (TYPE "2")? If the answer is "1", DISSPLA queries the user about terminal model

number, linespeed, resolution mode, etc. (The user should refer to the DISSPLA documentation [3] for the correct answers to these questions.) DISSPLA plots will now appear on the terminal instead of line

printer plots. If the answer is "2", line printer plots (the default) will

be displayed on the terminal, and DISSPLA plots will be saved in the

local file ZZZPLOT for later access.

4) ENTER LABEL FOR INDEPENDENT VARIABLE 1 (MAX.

OF 47 CHARACTERS): If the model has only one independent

variable, the user enters a label for that variable to be used in DISSPLA plots (up to 47 characters). If the problem has more than one

independent variable, the user enters the label for the first one, and an

opportunity to enter labels for the others will occur later (see prompt

26)).

5) ENTER LABEL FOR DEPENDENT VARIABLE (MAX. OF

47 CHARACTERS): The user enters a label for the dependent

variable to be used in DISSPLA plots (up to 47 characters).

6) DO YOU WANT LINE PRINTER OUTPUT? If the answer is "Y" , the program writes hardcopy output to a local file named LPOUT

which is assigned to unit 10. This output includes line printer graphs

of the plots requested, as well as iteration details, tables of results and

all other information requested by the user.

7) ARE YOUR DATA, WEIGHTS, INITIAL VALUES, ETC. ALREADY IN A VARPRO-DEFINED FILE? In any session the

user may save the data, weights and initial values in a file (see prompt

21)). If such a file has been written during a previous session, the

user may enter "Y" and that file will provide the data, weights and

initial values for this session. The user is asked to enter the name of

the file (maximum of 26 characters) and is given the option to have

these initial values printed. The program then advances to prompt

21). If no such file exists, the user should respond with a

tt N" and the

subroutine will continue with prompt 8).

43
8) NO. OF DATA POINTS: The user enters the number of observations

(maximum of 300).

9) NO. OF LINEAR PARAMETERS: The user enters the number of

linear parameters (maximum of 10).

10) NO. OF NONLINEAR PARAMETERS: The user enters the

number of nonlinear parameters (maximum of 8).

11) NO. OF INDEPENDENT VARIABLES: The user enters the

number of independent variables per observation (maximum of 6).

12) ARE DATA IN A FILE? The data can be read in from a file con- taining one line per observation. Each line should contain the values

of the independent variables (max. of 6) followed by the corresponding value of the dependent variable. If the user answers "Y" to this

prompt, the program asks for the name of the file (max. of 26 characters), and, upon receiving it, reads the data into the program. The

program then skips to prompt 14).

13a) ARE WEIGHTS TO BE ENTERED WITH DATA? This

prompt is issued only in the user responded negatively to the one

in prompt 12). Such a response means that he must type in the data

from the terminal. For each observation he must type in the values of

the independent and dependent variables. He is also given the option

of typing in the weight to be associated with each observation. If the answer to this prompt is "Y", weights will be entered by the user. 13b) DATA POINTS: LIST INDEPENDENT VARIABLES AND

DEPENDENT VARIABLE (AND WEIGHTS), SEPARATED

BY BLANKS OR COMMAS An index number is printed for each

point. The user types in all the necessary values for each observation.

Values are accepted in free format with or without decimal points or exponents (e.g., 100, 100.1 or 1.001E+02 are all valid).

14) ARE WEIGHTS TO BE CALCULATED IN SUBROUTINE

WT? If the answer is "Y", the user must have previously written the

subroutine WT to calculate the weights. If the answer is not "V, all weights are set to 1.0 (equal weighting).

15) DO DEPENDENT OR INDEPENDENT VARIABLES NEED

TO BE TRANSFORMED BEFORE ANALYSIS? If the an- swer is "Y", the user must have previously written the subroutines

TRANS and BAKTRN to do the desired transformation and backtransformation (cf. Chapter 3). (Please note: If data is transformed,

44
all subsequent tables and graphs will contain the transformed values,

unless specifically stated otherwise.)

16) INITIAL VALUES OF NONLINEAR PARAMETERS An in- dex number appears for each nonlinear parameter. The user enters the initial estimate of that parameter.

17) HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER j: This prompt is repeated for each of the

nonlinear parameters otj. A term depends on a nonlinear parameter

if that parameter is contained in the formula that defines the term.

The user enters the number of dependent terms. The program then

queries "WHICH TERMS ARE THEY?". The user enters the in- dex numbers of all the dependent terms for that parameter, separated

by blanks or commas.

18) ALL DERIVATIVES IN CLOSED FORM? If the answer is "Y",

the user must have previously written the subroutine DERIV to calculate the derivatives. The program then advances to prompt 21). If the answer is not "Y", numerical derivatives must be specified (prompts

19) and 20).

19) ALL NUMERICAL DERIVATIVE INCREMENTS THE SAME?

IF the answer is "Y", the user is asked to enter the increment value

and the program advances to prompt 21).

20) NUMERICAL DERIVATIVE INCREMENTS An index number

appears for each increment value. The increments are ordered by

nonlinear parameter and then by term of the model. For example,

if terms 1 and 2 depend on nolinear parameter 1, and terms 2 and

3 depend on nonlinear parameter 2, then the increments would be

ordered DPH(1,1), DPH(2,1), DPH(2,2) and finally DPH(3,2). (The

first subscript corresponds to term and the second to parameter.) Any

single derivative may still be calculated in closed form in subroutine

DERIV by setting that increment value equal to zero. 21) WOULD YOU LIKE TO SAVE THE DATA JUST INPUT?

If the answer is "Y", the user enters a file name (maximum of 26

characters) for a local file into which the data just specified are read.

After the INVAR session is completed, the user can store this file as a permanent file which can then be used in any subsequent session to

automatically read in the data, weights and initial values (cf. prompt

7)). If the user chose to transform the data (prompt 15)), then this

file will contain transformed, not original, values.

45
22) DOES DATA FOLLOW AN EVENLY SPACED TIMES SERIES? If the answer is "Y", all line printer graphs containing the

independent variable as the X axis will be printed as vertical plots.

DISSPLA graphs containing the independent variable as the X axis

will have a longer X axis than the other DISSPLA graphs in this case.

7.3 The Command Point

23) NEXT STEP (TYPE "H" FOR HELP) This is the command

point where the user must decide what is to be done next. A list of

all the options is given in Table 7.1. Entering an "H" at the command

point will produce the following message:

THE COMMAND CODES FALL INTO FOUR MAIN CATEGORIES

PRINT- TO EXAMINE CURRENT VALUES OF INITIALIZED DATA

CHANGE- TO CORRECT ANY INITIALIZED VALUE, TRANSFORM

DATA. AND EITHER ADD OR DELETE DATA POINTS

RUN- TO RUN THE LEAST SQUARES PROGRAM TO FIND THE

OPTIMAL PARAMETERS, OR TO RUN A STATISTICAL

ANALYSIS OF THE DATA RETURNED

SAVE/ TO SAVE CURRENT DATA OR

STOP- TO STOP THE PROGRAM, WITH OR WITHOUT

SAVING THE INITIALIZED DATA.

TYPE "P'V'C'V'R", OR "S" FOR THE SPECIFIC COMMANDS IN

THE ABOVE CATEGORIES. TYPE "A" TO OBTAIN A LISTING OF

ALL COMMAND CODES.

Bascially the options fall into four categories:

1. print any of the current values,

2. change any of the current values or transform the data;

3. run the least squares program or the statistical analysis; or

4. save the current values and/or end the terminal session.

Whichever command is chosen, the program cues the user to give

all the necessary information to successfully execute that command.

After execution is complete, the program returns to the command

46
Table 7.1: The INVAR command codes

COMMAND CODES FOR PRINTING

PM

PN

PKG

PIV

PTW

PALF

PNDI

PINC

LP

NP

TP

NTP

PRINT THE NO.

PRINT THE NO.

PRINT THE NO.

PRINT THE NO.

OF DATA POINTS

OF LINEAR PARAMETERS

OF NONLINEAR PARAMETERS

OF INDEPENDENT VARIABLES

PRINT THE DATA POINTS ft WEIGHTS

PRINT THE VALUES OF THE NONLINEAR PARAMETERS

PRINT THE NUMERICAL DERIVATIVE INCREMENTS

PRINT WHICH TERMS DEPEND ON NONLIN PARAMETERS

TURN ON THE LINE PRINTER FLAG

TURN OFF THE LINE PRINTER FLAG

TURN ON THE FLAG TO PLOT ON CRT

TURN OFF THE FLAG TO PLOT ON CRT

COMMAND CODES FOR CHANGING

CM CHANGE NO. OF DATA POINTS

CN CHANGE NO. OF LINEAR PARAMETERS

CKG CHANGE NO. OF NONLINEAR PARAMETERS

CIV CHANGE NO. OF INDEPENDENT VARIABLES

CT CHANGE ONE OR MORE DATA POINTS

CW CHANGE ONE OR MORE WEIGHTS

CALF CHANGE ONE OR MORE NONLINEAR PARAMETERS

CNDI CHANGE ONE OR MORE NUMER DERIV INCREMENTS

CTRN TRANSFORM THE DATA

CINC CHANGE WHICH TERMS DEPEND ON NONLIN PARAMETERS

COMMAND CODES FOR RUNNING

RV RUN THE LEAST SQUARES PROGRAM

RS RUN A STATISTICAL ANALYSIS OF RESULTS

COMMAND CODES TO SAVE DATA FILE AND/OR STOP

SC SAVE CURRENT VALUES OF DATA ft CONTINUE

SD SAVE CURRENT VALUES OF DATA ft STOP

SR SAVE CURRENT VALUES OF INDEP. VARIABLES AND

RESIDUALS ft CONTINUE

SO STOP WITHOUT SAVING

47
point (unless, of course, the user requests the commands SO or SD,

which result in program termination).

The print, change and save/stop commands are self explanatory and

easy to follow, and will not be explained in detail here. The run commands take the user through one of two paths, both terminating

back at the command point, prompt 23).

7.4 Running The Least Squares Program

The code for running the least squares code VARPRO is "RV. It activates

the following prompts:

24) DO YOU WANT TO CHANGE ANY DEFAULTS IN THE

LEAST SQUARES PROGRAM? A list of default values appears.

If the answer is "Y", a list of codes to change these defaults is printed.

(For an explanation of the output printing index, see IPRINT in Appendix C.) The user enters the desired code and the program then

specifies the limits to which that criterion can be changed and asks

the user for a new value. The user enters the new value and is then

given the option to change any of the other defaults. Finally, the least

squares program begins.

25) INDICATE WHICH OF THE FOLLOWING SHOULD BE

PRINTED BY TYPING "Y": After the least square calculations

are completed, a list of 3 possible tables to be printed appears. After

each table is listed, the user enters "Y" if that table is to be printed, or any other character if that table is to be skipped. All selected tables are then printed.

26) PLOT ROUTINE CURRENTLY SET NOT TO DO ANY

PLOTS. DO YOU WANT TO START PLOTTING? This

prompt will be repeated for each of the independent variables. If there are more than one independent variable, the user will be asked

to type in a label for each one chosen except for the first whose la- bel was typed in response to prompt 4). When no more plots are

desired, the user responds to this prompt with something other than

"Y", and the program returns to the command point, prompt 23). If the response to this prompt is "Y" , a list of available plots appears,

and the user enters the number of the plots to be printed. The plot

of observed and predicted values is printed twice when the data has

48
been transformed, once using the transformed data and once using

the back-transformed data. The user also has the option of printing a

table of the back-transformed values.

7.5 Running the Statistics Programs

When a VARPRO fit has been completed and the program has returned to the command point, the user can obtain a statistical analysis of the fit and

the residuals by responding to the "NEXT STEP?" prompt with "RS"

which activates the following prompt:

27) INDICATE WHICH OF THE FOLLOWING SHOULD BE

PRINTED: A list of tables, data and plots to be printed appears.

After each item is listed, the user enters "Y" to have that item printed,

or some other character to skip it. All the selected items are then

printed. A more detailed explanation of the various statistical quantities is given in Chapter 8. The plot of standardized residuals vs. independent variable is printed once for each independent variable.

After all items are printed, the program returns to the command point,

prompt 23).

7.6 The VARPRO Error Diagnostic Messages

The actual least squares calculations are done in the VARPRO subroutine

package which checks for various errors and indicates their occurence with one of the following diagnostic messages.

PROBLEM TERMINATED FOR EXCESSIVE ITERATIONS

The maximum number of iterations was exceeded before the nonlinear

parameters converged to the accuracy requested. The easiest thing to do about this message is simply to start the iteration again by typing in an

U RV" command. The iteration will begin again using the nonlinear

parameter values obtained in the previous run as the initial estimates.

If the problem persists, the user may want to relax the convergence

tolerance before running the problem again, or examine the residuals

and fitted function for the unconverged parameter values.

PROBLEM TERMINATED BECAUSE OF ILL-CONDITIONING

This message is usually obtained when the iteration has found an apparent minimum at a point where the parameters are highly correlated

49
with one another. This situation can be verified by examining the cor- relation matrix. The occurence of this situation means that many

other parameter combinations will give a fit that is essentially identical to the one obtained. The fit is usually good, but the user should

be aware of the possibility that he has specified a model with more

nonlinear parameters than the data can support. It is our impression

that this message is overly pessimistic and can often be ignored.

INPUT ERROR IN PARAMETER N,KG,M,NPP2, OR MMAX.

One or more of the indicated parameters is in error. The parameter

values should be:

N The number of linear parameters,

KG The number of nonlinear parameters,

M The number of observations (data points),

NPP2 N + P + 2, where P is the number of nonvanishing partial

derivatives dQ/da,

MMAX The maximum number of observations allowed by the code.

Note that MMAX is set by parameter statements in the program.

ERROR - INC MATRIX IMPROPERLY SPECIFIED OR DISAGREES WITH NPP2. The INC matrix is an "incidence matrix"

that tells the program which of the $ functions depends on which of

the nonlinear parameters a. The user can print out this matrix by

typing "PINC" and change it by typing "CINC".

ERROR - WEIGHT( mm) NOT NONNEGATIVE. The weight

Wnnn is negative.

CATASTROPHIC FAILURE - COLUMN mm IS ZERO. SEE

RALL AND FUNDERLIC'S DOCUMENTATION. VARPRO

works with an M x N matrix A whose columns are set to the values

of the $ functions and their partial derivatives d$>/da evaluated at the observation points. The ordering of the columns is:

1. The first N columns contain the #jt functions that have linear

multipliers C^, ordered in the same way as they are in user subroutine CALC. Recall that any 3> which do not depend on the

nonlinear parameters should come first. Such function are called

"constant functions" by the VARPRO comment cards.

50
2. The (N+l)st column contains the 3>at+i which does not have a

linear coefficient C. If the model does not have a &n+i term,

then the (N+l)st column is left blank.

3. The final P columns contain the nonvanishing partial derivatives

dQic/dftj, ordered first on the a*, i.e., all of the d<&k/doL\ come

first, then all of the dQk/d<X2, and so on.

This error is usually caused either by a programming error in CALC

or DERIV, or because the iteration has produced a set of nonlinear

parameters a which yield vanishingly small values for one of the $

functions or one of the partial derivatives. An example of the second

situation might be an exponential <&(<) = exp(at) that gets out of hand

because the iteration produces very negative values of a. In some cases

this kind of error can be avoided by starting again with different initial estimates for the parameters, but if the iteration persists in converging

to this error, the user should carefully examine his model and/or his

CALC and DERIV subroutines.

ERROR- CONSTANT COLUMN MUST BE COMPUTED ONLY

WHEN ISEL = 1. See the explanation of the preceding error stop

for the definition of a constant column. We have never observed this error stop. It is caused when a constant column is computed more

than once. Users who get this stop should check for an error in his

subroutine CALC or in the incidence matrix INC.

51
Chapter 8

STATISTICAL

CONSIDERATIONS

This chapter is reprinted with slight alterations from [12, pages 54-58].

A variety of statistics are available after the optimal parameter values have been

calculated. These will be discussed in the order that they appear in the program.

One should bear in mind that for nonlinear least squares many statistics can only

be approximated, and their values should be treated accordingly.

The first statistical option is a table of the independent variables, observed

and predicted values of the dependent variable, weighted residuals and weights.

The weighted residual is the result of multiplying the residual by the appropriate

weighting factor, i.e., Wi {yi -rli ) = Wi ri

Recall that WEIGHT(I) = W?.

The variance-covariance matrix provides an estimate of the variance (or covariance) of the calculated parameters. Thus the diagonal elements are the variances

of the parameters (i.e., the average of the squared deviations). The off-diagonal

elements are the covariances of the parameters, where the covariance of two pa- rameters i and y is defined as the average cross-product of the i and y deviations

from their means x and y [5]. The covariance matrix is calculated as described by

Lawson and Hanson [10], from the formula

V = a

2 T-*T- T = a2 {TTT)~ 1

,

where the superscript T is transpose and a

2

is a scalar factor such that

2 {norm of residual) 2 a =

[degrees of freedom)

52
[The quantity (degrees of freedom) is the number of observations minus the number of unknown parameters.) The matrix T is the final reduced form of the matrix

of functions $ and derivatives d$/da used in calculating the optimal parameters.

The correlation matrix is a scaled version of the covariance matrix. Its x, y element is cov(x, y)/(a(x)a(y)), where a(i) represents the standard deviation of

the tth parameter. Thus the diagonal elements are all ones and the off-diagonal

elements are the correlation coefficients and have absolute value between zero and

one. The closer the values are to one, the stronger the relationship is between

the two parameters. A negative coefficient indicates an inverse relationship, i.e., y decreases as x increases; a positive coefficient indicates that x and y tend to increase (or decrease) together. Plots of the residuals and of the observed and

predicted values versus the independent variable give a pictorial view of the data

in the table. One should note, however, that the plot is of the actual residuals, and

not of the weighted residuals. On execution of the RS command, several other statistical options are available. The first is a table of the parameters, their standard deviations and the t-ratio (the

parameters divided by their standard deviations). The t-ratio follows the Student's

t-distribution if the sample population is normal or near-normal [5], It can be

used to determine the level of significance of the hypothesis that the mean of the

population takes on a particular value (e.g., the parameter takes on a particular

value) or to construct a confidence interval for the mean of a normal population

(e.g., the parameter) [2]. The weighted residual sum of squares (RSS) is the quantity which is minimized,

i.e., the sum of the squares of the weighted residuals. The residual mean square

(RMS), sometimes called the estimated variance, is the residual sum of squares divided by the degrees of freedom. The standard error is the square root of the

residual mean square. The coefficient of determination (also referred to as R2 or the square of the

multiple correlation coefficient) is another measure used to determine the goodness

of fit. It is calculated from the formula

rf = l

( RSS)> (CTSS) '

where the corrected total sum of squares is found by

m

(ctss) = £ w?(yi - y)

2

,

where y is the weighted mean of the observed values

The coefficient of determination is thus a normalized residual sum of squares and

will always have a value between zero and one. The closer the value is to one, the

53
better the fit. For example, a coefficient of .995 would indicate that 99.5 percent

of the corrected total sum of squares is accounted for by the fitted equation [5]. One also has the option of printing a table of values of the independent variables,

the observed and predicted values of the dependent variable, the residual and the

standardized residual. The standardized residuals are the residuals divided by their standard deviations and are calculated from the formula

(SR)i =

l(RM'.

V *7

where «< is the standard deviation of the predicted value of the tth parameter [8]. This latter quantity is given by

Si = yJeJ*V*Tei

where V is the covariance matrix described above, # is the matrix whose columns

are the fitting functions $5 evaluated at the observation points, and €j is the vector with a one in the tth element and zeroes everywhere else. Thus s* is the tth diagonal

element of the matrix #V#T

Several plotting features are also available. The standardized residual is used

to provide a common scale and to avoid distortion caused by differences in the

standard deviations of the residuals. The first is a plot of the standardized residuals versus the independent vari- able (s). If the values do not appear to lie in a vertical band about the axis, one

of several things may be indicated. If the variance steadily increases (or decreases) over the range of values to form a wedge-shaped graph, the implication is that

a suitable transformation or a weighted least squares analysis should have been

performed. If the values lie in a linear band that is not vertical, a linear term

should have been included in the model. If the band of values is curved, linear and

quadratic terms should have been included in the model. One should realize that

combinations of these defects can also occur [4]. The standardized residuals are also plotted against the predicted values. Here

again the analysis is similar to that above. Nonrandomness might indicate the need

to use weights or the omission of some important term from the model equation

[8]- A probability plot of the standardized residuals versus the expected value of

the standardized residuals is also provided [8]. In making this plot, it is assumed

that the measurement errors are normally and independently distributed. If the

model adequately represents the data, the points will lie roughly on a diagonal line extending approximately between the points (—3, —3) and (3,3) and will tend to be clustered toward the center point (0,0).

54
Chapter 9

A TIME SERIES

EXAMPLE - THE

WILDLIFE CYCLE

9.1 The Model and the Data

The populations of many species of wildlife in the northern latitudes of the

Northern Hemisphere exhibit a mysterious but very regular cycle of length

approximately 10 years [6], [9]. Figure 9.1 is a plot of the annual total

number of Canadian lynx pelts purchased from trappers in the MacKenzie

River district by the Hudson's Bay Company in the years 1821-1934. The

individual yearly totals are connected by straight line segments to emphasize

the time-series nature of the data. The 10-year cycle is very pronounced, and

there is a hint that it was modulated with a cycle of period approximately

40 years. Since population processes tend to be exponential in nature, the usual

procedure is to work with the natural logarithm of the population measure.

Figure 9.2 was obtained by taking the natural logarithms of the numbers

plotted in Figure 9.1 and then by subtracting a least squares trend line

y{t) = d + C2

t .

The purpose of the detrending was to facilitate a spectral analysis of the

time series. Figure 9.3 is the periodogram of the data plotted in Figure 9.2.

It is dominated by the "10-year" cycle which whose length appears to be

55
MRCKENZIE RIVER LYNX CRTCH

& ^

•& i

YEAR

Figure 9.1: The MacKenzie River lynx catch in the years 1821-1934.

56
DETRENDED LOGRRITHM OF LYNX CRTCH

Figure 9.2: Natural logarithm of the MacKenzie River lynx catch with the

least squares trend line removed.

57
DETRENDED LOGRRITHM OF LYNX CRTCH

PERIODOGRRM

FREQUENCY

Figure 9.3: Periodogram of the detrended natural logarithm of the MacKenzie River lynx catch

58
closer to 9.5 years than to 10 years. There is also a hint of a cycle of length

approximately 40 years.

Figures 9.2 and 9.3 suggest the following model for the MacKenzie River

lynx catch:

y(t) =Ci + C2 t + C3 sin

2tt

10

{t + ew) + C4 sin

J 40

where t is the year, y(t) is the natural logarithm of the lynx catch in that

year, T\q and T40 are the two periods (to be determined by the fit), and 0\q

and #40 are the corresponding phase constants. This model, as it is written,

has four linear constants Ci, C2 , C3, and C4, and four nonlinear constants

^10? ^40) ^10? and ^40- It is possible to use the trigonometric identity for the

sine of the sum of two angles to rewrite this model with 6 linear constants

and only two nonlinear constants, but we will leave it in the form above for

this example.

Taking

«i = 2\o , <*2 = T40 , "3 = #10 5 Q4 = #40 5 we can write the model in the form

r){t) = Ci*!(0 + C2 *2 (<) + C3 $3 («) + C4$4 («) ,

where

*i(t) = 1.0

$3 (0 = sin

*4(0 sin

This is the model that was used for making the fits in this example. The

initial estimates that were used are

a? = 9.0 «S

a _ 5.0 a _

38.0

15.0

The data used in this example were the numbers plotted in Figure 9.1,

i.e., the number of lynx caught as a function of year. These data were written as two columns of numbers in a file called LYNXCAT. Each line of the

file contained a year number and the total lynx catch for that year. The

model, however, applied to the natural logarithms of the catch numbers.

59
Therefore it was necessary to transform the data in the user-supplied suboutine TRANS before making the fit. When making this transformation,

the independent variable t was also transfomed by moving the origin to the

midpoint of the record (i.e., by subtracting 1878.0 from all of the year val- ues). This is a standard practice which improves the conditioning of the

fit.

9.2 The User-Supplied Subroutines

The user-supplied subroutines required to make the fit are given in the fol- lowing paragraph. These subroutines were all written into a single file, called

LYCALC, in accordance with the directions given in Chapter 3. Subroutine

CALC is quite straightforward. We decided to let the program compute numerical approximations to the derivatives of the functions 3>t with respect

to the nonlinear parameters Qy, so subroutine DERIV is only a dummy subroutine. We used equal weighting in making the fit, so subroutine WT is also a dummy subroutine. Subroutine TRANS makes the transformations

described at the end of the preceding section, and subroutine BAKTRN

inverts the transformation on the dependent variable.

The subroutines used for this example are:

CC

cc

c

c

c

c

SUBROUTINE CALC ( K, NL, N. MMAX . IV. C, ALF , T. ETA )

PARAMETER ( MAXN1=11 )

DIMENSION C(N), ALF(NL). T(MMAX.IV)

COMMON /BLKPH/ PH(MAXNl)

DATA PI /3. 141592654/

PH(1) = 1. PH(2) = T(K,1)

PH(3) = SIN((2*PI/ALF(1))*(T(K.1)+ALF(3)))

PH(4) = SIN((2*PI/ALF(2))*(T(K.1)+ALF(4)))

60
ETA = C(1)*PH(1)+C(2)*PH(2)+C(3)*PH(3)+C(4)*PH(4)

RETURN

END

C

SUBROUTINE DERIV ( K, T. MMAX . IV. ALF . NL. NP1 ) C

C

C ****** DUMMY SUBROUTINE ******

C

C

PARAMETER ( MAXN1=11, MAXKG=8 )

DIMENSION ALF(NL).T(MMAX,IV)

COMMON /BLKDPH/ DPH(MAXN1 .MAXKG)

DO 10 J=1,NL

DO 10 1=1. NP1

DPH(I.J) = 0.0

10 CONTINUE

C

C AT THIS POINT INSERT THE STATEMENTS NEEDED TO SET THE

C REQUIRED NON-ZERO DERIVATIVES IN THE ARRAY DPH(I.J).

C

C DPH(I.J) = PARTIAL DERIVATIVE OF THE ITH FUNCTION PHI

C WITH RESPECT TO THE JTH NON-LINEAR PARAMETER ALF

C EVALUATED AT OBSERVATION NUMBER K.

C I.E. EVALUATED AT T(K.L). L=1.IV

C

RETURN

END

61
c

c SUBROUTINE WT ( M. WEIGHT. Y ) C

C

DIMENSION Y(M),WEIGHT(M)

C

C ****** DUMMY SUBROUTINE ******

C

C TO GET 1/Y**2 WEIGHTING. SIMPLY ACTIVATE THE NEXT 3 STATEMENTS

C

C DO 10 1=1,

C WEIGHT(I) = 1.0/Y(I)**2

C 10 CONTINUE

C

C OTHER WEIGHTINGS CAN BE GOTTEN BY CHANGING THE STATEMENT INSIDE

C THE LOOP.

C

RETURN

END

SUBROUTINE TRANS ( T, BT, Y, BY. IV, M, MMAX )

DIMENSION T(MMAX.IV). BT(MMAX.IV). Y(M) . BY(M)

DO 20 1=1.

DO 10 J=1.IV

BT(I.J) = T(I.J)

T(I.J) = T(I.J) - 1878.0

10 CONTINUE

BY(I) = Y(I)

Y(I) = ALOG(YU))

62
20 CONTINUE

RETURN

END

C

C

C

C

C

C

SUBROUTINE BAKTRN ( ETA, BETA, M )

DIMENSION ETA(M), BETA(M)

DO 10 1=1,

M

BETA(I) = EXP(ETAU))

10 CONTINUE

RETURN

END

These subroutines were written into a file named SNCALC which was compiled and loaded with the INVAR programs and then executed in a teminal

session to be described in the following section.

9.3 The Terminal Session

The terminal session output that follows was obtained by executing INVAR2

on the NBS CYBER 855, running from a Qume QVT-201 terminal in screen mode. Since the terminal was not a plotting terminal, the DISSPLA plots were done after execution on the QMS laser printer with address B225.

Before reading this example, the following points should be noted:

1. Since the terminal was running in screen mode, the INVOKE,INVAR

command generated a new screen shown between the horizontal lines at the beginning of the session. This was done by the interactive

NOS Procedure to quiz the user for the location of the user-supplied

63
subroutines, the address of the QMS printer for the plots and hardcopy

output, and the name of the data file.

2. Since the derivatives were to be computed numerically, it was necessary

to specify nonlinear parameter increment values to be used in the finite difference approximations. We chose to use the same increment value,

0.001, for all four parameters.

The VARPRO least squares iteration terminates with the warning,

"PROBLEM TERMINATED BECAUSE OF ILL-CONDITIONING",

but this does not indicate a fatal error. It occurs fairly often when

using VARPRO, but the linear and nonlinear parameters calculated

by the program are usually satisfactory.

Since this is a time-series problem, all the line-printer graphs with

time as the independent variable are plotted in a vertical rather than

a horizontal format.

After each line printer graph, the program stops with a question mark

(?) prompt to which we replied with a

U C followed by a carriage

return. The prompt was apparently issued by the TK4014 program

used at NBS to generate the DISSPLA plots for the laser printer.

Cyber users should always respond to this prompt with something more than just a carriage return in order to avoid an end-of-file error

stop.

The terminal session follows:

/INVOKE, INVAR

64
INTERACTIVE NONLINEAR LEAST SQUARES PROGRAM

CALCF FORTRAN file of user-supplied subroutines: LYCALC_

QMSADDR QMS Address for output (Default = NONE) : B225 DATAF Data file (Default = NONE) : LYNXCAT

INV102 INVAR1 or INVAR2 (Default = INVAR2)

Specify values and press NEXT when ready

INVARX

PROC . INVARX , CALCF [ , qMSADDR] [ , DATAF] [ , INV102]

Executes one of two versions, INVAR1 or INVAR2, of INVAR, an

interactive, variable separable, nonlinear least squares program.

Local file CALCF contains FORTRAN coded, user-supplied subroutines.

Hardcopy output and plots can be sent to QMS laser printer with

address QMSADDR.

Permanent file DATAF contains data to be fit (not needed if data

already in a local file, or is to be typed in during execution).

F3 HELP FWD F5 HELP F6 QUIT

ENTER TITLE (MAXIMUM OF 40 CHARACTERS)

? MACKENZIE RIVER LYNX CATCH

DO YOU WANT DISSPLA PLOTS? (TYPE "Y" FOR YES)

? Y

DO YOU WANT THESE PLOTS SHOWN INTERACTIVELY ON A PLOTTING SCREEN (TYPE "1"); OR SAVED IN A FILE TO BE PRINTED OUT LATER (TYPE "2")?

? 2

ENTER LABEL FOR INDEPENDENT VARIABLE (MAX. OF 47 CHARACTERS)

? YEAR

ENTER LABEL FOR DEPENDENT VARIABLE (MAX. OF 47 CHARACTERS):

? CATCH

DO YOU WANT LINE PRINTER OUTPUT?

? Y

ARE YOUR DATA, WEIGHTS. INITIAL VALUES, ETC. ALREADY IN A VARPRO-DEFINED FILE ?

65
? N

NO. OF DATA POINTS:

? 114

NO. OF LINEAR PARAMETERS:

? 4

NO. OF NONLINEAR PARAMETERS:

? 4

NO . OF INDEPENDENT VARIABLES

:

? 1

ARE DATA IN A FILE?

? Y

WHAT IS THE FILENAME?- -NO MORE THAN 26 CHARS

(TYPE "QUIT" TO EXIT)

? LYNXCAT

ARE WEIGHTS TO BE CALCULATED IN SUBROUTINE WT?

? N

WEIGHTS ARE SET TO 1

DO DEPENDENT OR INDEPENDENT VARIABLES

NEED TO BE TRANSFORMED BEFORE ANALYSIS?

? Y

INITIAL VALUES OF NONLINEAR PARAMETERS

1 ? 9.0

2

? 38.0

3

? 6.0

4

? 15.0

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 1 ? 1

66
WHICH TERMS ARE THEY?

? 3

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 2:

? 1

WHICH TERMS ARE THEY?

? 4

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 3:

? 1

WHICH TERMS ARE THEY?

? 3

HOW MANY TERMS OF THE MODEL DEPEND ON NONLINEAR PARAMETER 4:

? 1

WHICH TERMS ARE THEY?

? 4

ALL DERIVATIVES IN CLOSED FORM?

? N

ALL NUMERICAL DERIVATIVE INCREMENTS THE SAME?

? Y

INCREMENT VALUES:

? .001

WOULD YOU LIKE TO SAVE THE DATA JUST INPUT?

? N

DOES DATA FOLLOW AN EVENLY SPACED TIME SERIES?

? Y

NOTE: THE COMMAND "SC" WILL SAVE YOUR DATA FILE AT ANY TIME

NEXT STEP? (TYPE "H" FOR HELP)

? RV

DO YOU WISH TO CHANGE ANY DEFAULTS IN THE LEAST SQUARES PROGRAM?

(50: MAI. NO. OF ITERATIONS)

(l.OE-6: CONVERG. TOLERANCE)

(1: OUTPUT PRINTING INDEX)

67
? N

1 MACKENZIE RIVER LYNX CATCH

WEIGHTED NORM OF RESIDUAL = 1 .3034692E+01

NU = .lOOOOOOE+01

ITERATION 1 NONLINEAR PARAMETERS

9.4014616E+00 3.7990261E+01

6.2013727E+00 1 .9306338E+01

1 WEIGHTED NORM OF RESIDUAL = 9.3597487E+00

NU = .BOOOOOOE+OO

NORM(DELTA-ALF) / NORM(ALF) = 9.862E-02

ITERATION 2 NONLINEAR PARAMETERS

9.6511209E+00 3 . 9674402E+01

4 . 9024596E+00 2 . 1067477E+01

1 WEIGHTED NORM OF RESIDUAL = 7 . 5367307E+00

NU = .2BOOOOOE+00

NORM(DELTA-ALF) / NORM(ALF) = 6.310E-02

ITERATION 3 NONLINEAR PARAMETERS

9.6094909E+00 4.0501357E+01

4.8074688E+00 2. 196903BE+01

1 WEIGHTED NORM OF RESIDUAL = 7 . 2985779E+00

NU = . 1250000E+00

NORM(DELTA-ALF) / NORM(ALF) = 2.597E-02

ITERATION 4 NONLINEAR PARAMETERS

9.6189932E+00 4.0636899E+01

4.792B146E+00 2. 1908877E+01

68
1 WEIGHTED NORM OF RESIDUAL = 7 . 2930963E+00

NU = .6250000E-01

NORM(DELTA-ALF) / NORM(ALF) = 1.362E-03

ITERATION B NONLINEAR PARAMETERS

9.6198493E+00 4.0566491E+01

4.7910923E+00 2 . 1917641E+01

1 WEIGHTED NORM OF RESIDUAL = 7. 2930368E+00

NU = .3126000E-01

NORM (DELTA- ALF) / NORM(ALF) = 6.B28E-04

ITERATION 38 NONLINEAR PARAMETERS

9.6199005E+00 4 .0579944E+01

4.7910060E+00 2. 1914852E+01

28 WEIGHTED NORM OF RESIDUAL = 7 . 2930340E+00

NU = B548352E+02

NORM(DELTA-ALF) / NORM(ALF) = 7.510E-10

ITERATION 39 NONLINEAR PARAMETERS

9.6199005E+00 4 .0579944E+01

4.7910060E+00 2 . 1914862E+01

1 WEIGHTED NORM OF RESIDUAL = 7 . 2930340E+00

STEP RETRACTED. NU = 8.322B286E+01

ITERATION 40 NONLINEAR PARAMETERS

9.6199006E+00 4.0679944E+01

69
4.7910060E+00 2.19148B2E+01

2 WEIGHTED NORM OF RESIDUAL = 7 . 2930340E+00

STEP RETRACTED, NU = 1 . 2483793E+02

PROBLEM TERMINATED BECAUSE OF ILL-CONDITIONING

LINEAR PARAMETERS

.6692901E+01 3384831E-02

NONLINEAR PARAMETERS

.9619901E+01 . 4057994E+02

WEIGHTED NORM OF RESIDUAL =

WEIGHTED ESTIMATED VARIANCE =

. 1436036E+01 .4943990E+00

4791006E+01 .2191486E+02

.7293034E+01

.B017768E+00

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED

BY TYPING "Y": TABLE OF OBSERVED ft PREDICTED VALUES, WTD RESIDUAL:

? Y

COVARIANCE MATRIX:

? N

( ? Y

1

CORRELATION MATRIX:

DEX IND VAR 1 OBS VAL PRED VAL WTD RESID WEIGHTS

1 -5 70000E+01 5 59471E+00 6 2373BE+00 -6 42634E-01 1 OOOOOE+OO

2 -6 60000E+01 5 77144E+00 6 63420E+00 1 37242E-01 1 OOOOOE+OO

3 -5 50000E+01 6 37161E+00 5 BB072E+00 8 20896E-01 1 OOOOOE+OO

4 -5 40000E+01 6 76964E+00 6 03666E+00 7 33093E-01 1 OOOOOE+OO

6 -5 30000E+01 7 29641E+00 6 90279E+00 3 93622E-01 1 OOOOOE+OO

6 -6 2OOO0E+O1 7 9448BE+00 7 799B0E+00 1 4B343E-01 1 OOOOOE+OO

7 -6 10000E+01 8 27B89E+00 8 3B964E+00 -8 37B38E-02 1 OOOOOE+OO

8 -6 OOOOOE+01 8 68997E+00 8 3B012E+00 3 3984BE-01 1 . OOOOOE+OO

70
107 4.90000E+01

108 6.00000E+01

109 6 . 10000E+01

110 B.20000E+01

111 6.30000E+01

112 6.40Q00E+01

113 5.B0000E+01

114 B.60000E+01

7 . 33759E+00

6 . 27099E+00

6.18415E+00

6.49627E+00

6 . 90776E+00

7.37149E+00

7.88495E+00

8.1303BE+00

6.68261E+00

B.01960E+00

6.02B34E+00

B.60B13E+00

6.63246E+00

7.44190E+00

7.97913E+00

7.94642E+00

1.7BB08E+00

1.2B139E+00

1.1B881E+00

8.90136E-01

3.7B292E-01

-7.04083E-02

-9.4181BE-02

1 . 83936E-01

1 . 0000OE+O0

1 . 00000E+00

1 . 00000E+00

1 . 00000E+00

1 . 00000E+00

1 . 00000E+00

1 . OO00OE+O0

1 . 00000E+00

CORRELATION MATRIX

ROW COLUMN COLUMN COLUMN COLUMN

1 1 OOOOOOE+00 6 913627E-03 9 202B10E-03 1 B8B469E-02

2 6 913627E-03 1 OOOOOOE+OO -7 093361E-02 2 470090E-01

3 9 202B10E-03 -7 093361E-02 1 OOOOOOE+OO 1 993396E-03

4 1 B8B469E-02 2 470090E-01 1 993396E-03 1 OOOOOOE+OO

6 2 1B9597E-02 4 740477E-02 -1 747922E-02 -1 876114E-02

6 4 382337E-02 -2 098704E-01 9 0258BBE-02 -8 623688E-02

7 -1 391344E-02 -2 037B83E-02 -4 B4B04BE-03 -B .482076E-03

8 1 001278E-01 -1 B43260E-01 B 609647E-02 -4 .3689B0E-02

ROW COLUMN COLUMN COLUMN COLUMN

1 2 1B9B97E-02 4 382337E-02 -1 391344E-02 1 001278E-01

2 4 740477E-02 -2 098704E-01 -2 037B83E-02 -1 B43260E-01

3 -1 747922E-02 9 02B8B6E-02 -4 B4S04BE-03 6 609647E-02

4 -1 87B114E-02 -8 623688E-02 -B 482076E-03 -4 3689B0E-02

B 1 OOOOOOE+OO -4 964643E-02 1 177933E-01 -1 .404177E-03

6 -4 964643E-02 1 OOOOOOE+OO 3 26B023E-02 B 878888E-01

7 1 177933E-01 3 266023E-02 1 OOOOOOE+OO 9 272713E-03

8 -1 .404177E-03 B .878888E-01 9 272713E-03 1 OOOOOOE+OO

PLOT ROUTINE CURRENTLY SET NOT TO DO ANY PLOTS.

DO YOU WANT TO START PLOTTING?

? Y

IF YOU WANT TO PLOT RESIDUALS, TYPE "1"

IF YOU WANT TO PLOT OBSERVED * PREDICTED VALUES, TYPE "2' IF YOU WANT TO PLOT BOTH. TYPE "3":

? 3

71
L RESIDUALS

IND VAR 1

MEAN IS .47144386E-10

GRAPH INTERVAL IS .6748B643E-01

-COMPUTED FROM MAXIMUM AND MINIMUM VALUES

-1.62 -.94 -.27 .41 1.08 1.76

-1.28 -.61 .07 .74 1.42

j 1 J j J 1 1 x j j.

I

1 I * 2 I 3 I 4 I

5 I 6 I 7 I 8 I 9 I 10 I 11 I 12 I 13 I 14 I 15 I 16 I 17 I 18 I 19 I 20 I

- . 64263E+00

13724E+00

82090E+00

. 73309E+00

. 39362E+00

. 14534E+00

- . 83764E-01

. 33984E+00

. 73926E+00

. 10130E+01

. 32294E+00

- . 82240E+00

-.234B1E+00

-.38987E+00

- . 8468BE+00

. 13947E+00

. 24994E-02

. 52962E+00

68781E+00

. 17464E+00

95 I 96 I 97 I 98 I 99 I 100 I 101 I 102 I 103 I 104 I 105 I. 15581E+00

88637E+00

. 15373E+00

. 11610E+01

. 86253E+00

.84513E+00

. 82754E+00

. 11315E+01

. 69832E+00

. 24565E-01

. 83975E+00

72
106 I 107 I 108 I 109 I 110 I 111 I 112 I 113 I 114 I

I —I —I —I- ! 1 J J !-

.1B063E+01

.17B51E+01

.12B14E+01

.11B88E+01

. 89014E+00

.37529E+00

-.70408E-01

-.94181E-01

. 18394E+00

? C

I OBSERVED (*) ft PREDICTED (+) VALUES (TRANSFORMED DATA)

IND VAR 1

XI VECTOR

MEAN IS 6.6798163

X2 VECTOR

MEAN IS 6.6798163

GRAPH INTERVAL IS .10377634

-COMPUTED FROM MAXIMUM AND MINIMUM VALUES

3.66 4.70 B.74 6.78 7.81 8.8B

4.18

I 1 1-

I

5.22 6.26

--I 1— -I7 . 30 8 . 33

XI = * X2 = +

1 I * + . 559E+01 .624E+01

2 I + * .B77E+01 .563E+01

3 I + * 637E+01 .BB5E+01

4 I + * 677E+01 604E+01

5 I + * [ .730E+01 .690E+01

6 I + * [ .794E+01 .780E+01

7 I * + [ . 828E+01 .836E+01

8 I + * [ .869E+01 .83BE+01

9 I + * C .8B1E+01 .777E+01

10 I + * [ .78BE+01 .684E+01

11 I + * C .626E+01 .B94E+01

12 I * + C .4B8E+01 .B41E+01

13 I * + [ .B21E+01 .B46E+01

14 I * + [ .563E+01 .602E+01

15 I * + t .601E+01 .686E+01

16 I +* ; [ .773E+01 .759E+01

17 I 2 [ .790E+01 .789E+01

18 I + * [ .813E+01 .760E+01

19 I + * [ .751E+01 .682E+01

20 I + * [ .601E+01 .B84E+01

21 I 2 [ .502E+01 .B04E+01

73
22 I * 23 I 24 I 25 I 26 I 27 I

381E+01

422E+01

536E+01

.630E+01

.694E+01

766E+01

.473E+01

.602E+01

.678E+01

.668E+01

.736E+01

.752E+01

93 I 2 ] .824E+01 .818E+01

94 I * + [ . 804E+01 .847E+01

96 I * + [ . 800E+01 .816E+01

96 I + * [ .824E+01 .735E+01

97 I + * [ .651E+01 . 636E+01

98 I * + [ .439E+01 . 656E+01

99 I * + [ .438E+01 .624E+01

100 I * + [ .468E+01 .563E+01

101 I * + [ . 643E+01 .626E+01

102 I * + t . 599E+01 .712E+01

103 I * + [ .703E+01 .773E+01

104 I 2 [ .780E+01 .782E+01

106 I + * [ 818E+01 .734E+01

106 I + * [ .798E+01 .648E+01

107 I + * [ .734E+01 .558E+01

108 I + * [ .627E+01 602E+01

109 I + * [ 618E+01 . 503E+01

110 I + * C .650E+01 .561E+01

111 I + * [ .691E+01 .653E+01

112 I 2 [ .737E+01 .744E+01

113 I * + [ .788E+01 .798E+01

114 I + v ;[ .813E+01 .796E+01

I-- __! ! 1 J.. —I—-I 1- —I— --i —

? c

OBSERVED (*) ft PREDICTED ( + ) VALUES (BACK-TRANSFORMED DATA) IND VAR 1

XI VECTOR

MEAN IS 1522.7106

X2 VECTOR

MEAN IS 1346.4738

GRAPH INTERVAL IS 139.04000

-COMPUTED FROM MAXIMUM AND MINIMUM VALUES

39.00 1429.40 2819.80 4210.20 6600.60 6991.00

74
I I

734

• 20 2124.60 3516.00 4905.40 6295.80

.t j j J J J 1 J j ] XI = * 12 = +

1 I * + . 269E+03 .611E+03

2 I 2 .321E+03 . 280E+03

3 I + * .685E+03 .257E+03

4 I + * .871E+03 .418E+03

5 I + * . 148E+04 . 996E+03

6 I + * . 282E+04 . 244E+04

7 I * + [ .393E+04 . 427E+04

8 I + * C . 694E+04 .423E+04

9 I + * t .495E+04 . 236E+04

10 I + * [ .268E+04 . 936E+03

11 I + * [ . 623E+03 .379E+03

12 I* + [ .980E+02 . 223E+03

13 I 2 [ . 184E+03 . 233E+03

14 I * + [ . 279E+03 .412E+03

15 I * + [ . 409E+03 .954E+03

16 I + * [ .229E+04 . 199E+04

17 I 2 ][ . 269E+04 .268E+04

18 I + * [ .341E+04 .201E+04

19 I + * I[ .182E+04 .917E+03

20 I + * [ . 409E+03 . 343E+03

21 I 2 [ .151E+03 .154E+03

22 I* + [ .45OE+02 .113E+03

23 I* + [ .680E+02 .152E+03

24 I * + [ . 213E+03 .324E+03

25 I * + [ . 546E+03 . 800E+03

26 I * + [ . 103E+04 . 156E+04

27 I + * [ .213E+04 . 184E+04

28 I + * [ .254E+04 .122E+04

29 I + * [ .957E+03 . 547E+03

30 I + * [ .361E+03 . 232E+03

90 I + 91 I 92 I 93 I 94 I 95 I 96 I 97 I 98 I* + 99 I*+

+ *

* +

* +

808E+03

139E+04

. 271E+04

380E+04

. 309E+04

. 299E+04

379E+04

.674E+03

.810E+02

.800E+02

. 386E+03

.716E+03

. 170E+04

357E+04

.475E+04

. 349E+04

. 156E+04

. 578E+03

. 269E+03

. 190E+03

75
100 I* + 101 I * + 102 I * 103 I 104 I 105 I 106 I + 107 I + 108 I + * 109 I + * 110 I + * 111 I 112 I 113 I 114 I

j 1 1 j J ! J J j j__.

. 108E+03

.229E+03

.399E+03

.113E+04

. 243E+04

. 367E+04

. 294E+04

. 154E+04

. 629E+03

.48BE+03

.662E+03

. 100E+04

. 169E+04

. 266E+04

.340E+04

? C

DO YOU WISH TO PRINT A TABLE OF THE BACK- TRANSFORMED VALUES?

? Y

OBSERVED AND PREDICTED VALUES BACK-TRANSFORMED

DEX IND VAR 1 OBS VAL PRED VAL

1 1 82100E+03 2 69000E+02 6 11499E+02

2 1 82200E+03 3 21000E+02 2 79835E+02

3 1 82300E+03 5 8B000E+02 2 57422E+02

4 1 82400E+03 8 71000E+02 4 18447E+02

6 1 82600E+03 1 47600E+03 9 95049E+02

6 1 . 82600E+03 2 82100E+03 2 43939E+03

7 1 .82700E+03 3 92800E+03 4 27116E+03

8 1 82800E+03 5 94300E+03 4 23071E+03

107 1 92700E+03 1 53700E+03 2 6B738E+02

108 1 92800E+03 5 29000E+02 1 B13B1E+02

109 1 92900E+03 4 8B000E+02 1 B2222E+02

110 1 93000E+03 6 62000E+02 2 71817E+02

HI 1 93100E+03 1 00000E+03 6 87089E+02

112 1 93200E+03 1 B9000E+03 1 70B98E+03

113 1 . 93300E+03 2 6B700E+03 2 91940E+03

114 1 . 93400E+03 3 39600E+03 2 82643E+03

.2B1E+03

. B24E+03

. 124E+04

. 228E+04

. 249E+04

. 1B4E+04

.661E+03

. 266E+03

. 151E+03

.1B2E+03

. 272E+03

. 687E+03

. 171E+04

. 292E+04

. 283E+04

NEXT STEP? (TYPE "H" FOR HELP)

76
? RS

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED: PARAMETERS. THEIR STD DEV * T-RATIO:

? Y

RES SUN OF SQS, RES MEAN SQ, RES STD DEV * COEFF OF DETERMINATION:

? Y

TABLE OF OBS. PRED VAL. RESID ft STD RESID:

? Y

PLOT OF STD RES VS IND VAR:

? Y

PLOT OF STD RES VS PRED VAL:

? Y

NORMAL PROBABILITY PLOT:

? Y

PARAMETER S.D. OF PARM. T-RATIO

6.692901E+00

3.384831E-03

1 . 436036E+00

4.943990E-01

9.619901E+00

4.067994E+01

4.791006E+00

2.191485E+01

6.673607E-02

2.129423E-03

9.396270E-02

9.466946E-02

2.974943E-02

1.679599E+00

1.013214E-01

1.663618E+00

1.002891E+02

1.B895B3E+00

1.628467E+01

6 . 222370E+00

3.233642E+02

2.4160B0E+01

4.728623E+01

1.401548E+01

WTD RESIDUAL SUM OF SQUARES: 6.318834E+01

WTD RESIDUAL MEAN SQUARE: 6.017768E-01

WTD RESIDUAL STANDARD ERROR: 7.083621E-01

COEFFICIENT OF DETERMINATION ( R- SQUARE) : 7.128861E-01

INDEX IND VAR 1 OBS VAL PRED VAL RESID STD RES

77
1 -6 70000E+01 5 69471E+00 6 2373BE+00 -6 42634E-01 -9 3B023E-01

2 -6 60000E+01 B 77144E+00 6 63420E+00 1 37242E-01 2 04622E-01

3 -6 BOOOOE+01 6 37161E+00 5 BB072E+00 8 20896E-01 1 19547E+00

4 -6 40000E+01 6 76964E+00 6 0365BE+00 7 33093E-01 1 06327E+00

6 -5 30000E+01 7 29641E+00 6 90279E+00 3 93622E-01 6 88660E-01

6 -5 20000E+01 7 94485E+00 7 799B0E+00 1 4B343E-01 2 1BB03E-01

7 -B 10000E+01 8 27689E+00 8 36964E+00 -8 37638E-02 -1 22281E-01

8 -B OOOOOE+01 8 68997E+00 8 3B012E+00 3 39846E-01 4 96178E-01

107 4 90000E+01 7 337B9E+00 B B82B1E+00 1 7BB08E+00 2 61812E+00

108 6 00000E+01 6 27099E+00 6 01960E+00 1 2B139E+00 1 84149E+00

109 B 10000E+01 6 18416E+00 5 02B34E+00 1 1B881E+00 1 70112E+00

110 5 20000E+01 6 49B27E+00 B 60B13E+00 8 9013BE-01 1 323B4E+00

111 5 30000E+01 6 90776E+00 6 63246E+00 3 7B292E-01 5 63784E-01

112 5 40000E+01 7 37149E+00 7 44190E+00 -7 04083E-02 -1 0B332E-01

113 5 B0000E+01 7 8849BE+00 7 97913E+00 -9 4181BE-02 -1 40B90E-01

114 5 60000E+01 8 1303BE+00 7 94642E+00 1 83936E-01 2 80809E-01

ENTER

? 1 ANY INTEGER TO CONTINUE

STANDARDIZED RESIDUAL VS INDEPENDENT VARIABLE 1

MEAN IS .16607042E-02

GRAPH INTERVAL IS .99228060E-01

-COMPUTED FROM MAXIMUM AND MINIMUM VALUES

-2.34 -1.36 -.36 .63 1.63 2.62

-1.8B -.86 .14 1.13 2.12

j j J j j x x ! ! j...

I

1 I * 2 I 3 I 4 I 6 I 6 I 7 I 8 I 9 I 10 I

11 I 12 I

-.93502E+00

. 20462E+00

.11956E+01

. 10633E+01

. B8866E+00

.21650E+00

- . 12228E+00

.49518E+00

. 10903E+01

. 16080E+01

. 47927E+00

-.12156E+01

78
13 I 14 I IB I 16 I 17 I 18 I 19 I 20 I

-.34910E+00

-.68756E+00

-.12772E+01

.20741E+00

.36642E-02

.77381E+00

.10081E+01

. 26690E+00

95 I 96 I 97 I 98 I 99 I 100 I 101 I 102 I 103 I 104 I 10E I 106 I 107 I 108 I 109 I 110 I 111 I 112 I 113 I 114 I

? c

.j ! ! j j j j j j.

- . 22937E+00

.1316BE+01

.22844E+00

- . 17096E+01

-.12623E+01

-.12435E+01

-.12278E+01

- . 16747E+01

-.1022BE+01

-.35890E-01

.12423E+01

.22550E+01

.26181E+01

.18416E+01

.17011E+01

.13236E+01

. B6378E+00

- . 10B33E+00

- . 14059E+00

.28081E+00

STANDARDIZED RESIDUAL VS. PREDICTED VALUES

1 Y DIV= .99228E-01 AND 1 X DIV= .77613E-01

_! ! j ! ! j_

2.6181 - *

79
1.6258

.6336

.3687

1 3510

2*

i * * *

* * *

* ** * * * *

* * * * 2

* * * *** ** * * * * * * *

*

* 2 * *

* * ** 2 *

* * *

* *

* *

* * *

* * * *

* * * ** * * 2

* *

*

*

2.3433 -

I _! j J 4.7307 5.5069 6.2830

._! ! j_

7.0591 7.8363 8.6114

? C

80
NORMAL PROBABILITY PLOT

1 Y DIV= .99228E-01 AND 1 X DIV

1- I—

10128

--I-- 1—

2.6181

1-

* -

1.6268

.6336

.3587

1.3610

* 232

**

* 3* 3

24*

*23 64*

* 23

43

2 *2 22

*2

2

*2 *232*

* *2
*

* -2.3433 - *

I

_ x j ! j j j_

-2.6166 -1.6038 -.4909 .6219 1.6348 2.6476

? C

NEXT STEP? (TYPE "H" FOR HELP)

? SO

END OF DISSPLA 9.0 -- 10773 VECTORS IN 6 PLOTS. RUN ON 4/30/87 USING SERIAL NUMBER AT NBS/CS**2 - GAITHERSBERG. MD. PROPRIETARY SOFTWARE PRODUCT OF ISSCO, SAN DIEGO, CA. 3644 VIRTUAL STORAGE REFERENCES; 4 READS; WRITES.

9.4 The Graphical Output

The 6 DISSPLA plots generated by the run are shown in Figures 9.4-9.9. The last three of these were generated by the statistical analysis initiated

by the a RS" command. They provide valuable diagnostic information for

assessing the quality of the fit. The standardized residuals are the residuals divided by their standard deviations. The fact that no wedge-shaped

patterns appear in Figures 9.7 and 9.8 indicates that we were probably

right to use equal weighting in the VARPRO least squares procedure. The

distribution of the points in Figure 9.9 along a straight line through the

points (-3.0,-3.0) and (3.0,3.0) is also a good indicator, but it is obvious

from inspecting the residual time series (Figure 9.4) that the residuals are

significantly correlated. A really good model would give white-noise residuals.

82
en

_j

en

CD

CD

LJ

en

so/,

SldDQIS3d

Figure 9.4: The residuals for the fit (transformed data - model)

83
C_)

cr

X

en

uJ

en

N

LJ

Q

CDU

21

CCo

Li- en

CC

Figure 9.5: The transformed data (boxes) and the fitted model (continuous

line)

84
CJ

CE

CJ

X

~ZL

LJ

LJ

N

LJ

cr

CEQ

D

LJr

QCO

U_LOZ

CEITI—

i *:CECD

<k <?5 <?3 ^ 'OfS> % *Jr MJlov

Figure 9.6: The original data and the back-transformed fit

85
enenQ

en

Q

en

LJ

<r. (9. <?. siunais3y Q3ZiayyjNbis

Figure 9.7: Standardized residuals versus the independent variable

86
STRND. RESIDURLS VS. PRED. VRLUES

CATCH

Figure 9.8: Standardized residuals versus the model prediction

87
NORMRL PROBRBILITY PLOT

,v ,v <v v v

EXPECTED STANDARDIZED RESIDUALS

Figure 9.9: Standardized residuals versus the expected standardized residuals when it is assumed that the measurement errors are independently,

normally distributed

88
Chapter 10

AN EXAMPLE WITH

TWO INDEPENDENT

VARIABLES - VISCOSITY

OF ORGANIC

COMPOUNDS

10.1 The Model and the Data

As a final example we present a fit of a bivariate model to measured viscosity

data for cyclic polystyrene molecules [11]. The model is

V(T,M',a,b,a,bh T00 ) = (aM+6MQ )10^fe ,

where the measured temperature T (degC) and the molecular weight M

of the molecule are the indepedent variables of the fit, and the measured

viscosity V is the dependent variable. The constants a, 6, a, 6i, and T^ are the parameters to be determined by the fit. Since the viscosity data ranges

over 12 orders of magnitude, the fitting is much easier is we take log10 V as the independent variable for the fit. Thus the model becomes

r}(T,M) = \og l0 V = b1¥-y- + log 10 (aM + bMa

) ,

89
which can easily be written in the variable separable form described in Chapter 2. The model has a single linear parameter

C\ — b\ ,

and 4 nonlinear paramenters

ai = T<» , «2 = a , 0:3 = 6 , a.\ — a .

Thus the model can be written in the form

r,{T,M) = d*i(T, J^) + *2 {T,M) ,

where

&2 {T,M) - log 10 (a2 M + a3 Ma*) = log 10 (aM + bMa

) .

The data consisted of 159 observed triplets (Ti,Mj,Vi), contained in a

file named VISDAT. To make the desired fit, it was necessary to provide a subroutine TRANS to convert these observations to triplets of the form

(Ti, M{, r)i) and a subroutine BAKTRN to backtransform the fit once it had

been obtained. These subroutines are given in the following section. On

the first run of the program, the transformed data were saved, by using

the SD command (see Table 7.1), in a file called VSTRDAT, and this file was used as input for subsequent runs, including the one in the terminal

session to be given below. The input file contained the final estimates of

the nonlinear parameters, saved from the previous run. These estimates,

except for the one of 0:3 — 6, were used as initial estimates in the current

run. The statistical analysis from the previous run indicated a high degree

of uncertainty for the value of 03. The current run checks the stability of

the fit by changing the estimate of 03 and repeating the calculation.

10.2 The User-Supplied Subroutines

The user supplied subroutines required to make the fit are given in the

following paragraph. These subroutines were all written into a single file, called VSCALC, in accordance with the instructions given in Chapter 3. Subroutine CALC is a straightforward coding of the model equations given

90
in the preceding section. Subroutine DERIV calculates the derivatives with

respect to the parameters from the exact, closed form expressions:

d$! _ l a«i _

(T-Too)2 '

9q2 da$ da^ o,

!** = o,

da2 ~ lU*>10\c ) aM+bM<

d$i - ] na f r

\ Ma 9a 3 — lo&10V c y aM+6M'

9$2 9a 4 logl0le J aM+bM°

Subroutine WT is designed to give \/yf weighting, but it was not used

in the present calculations. The equal weighting default option was used

instead. If the user is not sure beforehand about which weighting should

be used, it is a good idea to include the WT subroutine given here in case

it should be needed. The interactive driver program allows the user to choose which weighting scheme to use at execution time or even to change

the weighting scheme if the first one chosen does not work satisfactorily.

Subroutines TRANS and BAKTRN make the transformations described in the preceding section.

The subroutines used for this example are:

C

C

C

C

C

C

CC

SUBROUTINE CALC ( K, NL, N, MMAX , IV. C, ALF , T, ETA )

PARAMETER ( MAXN1 = 11 )

DIMENSION C(N), ALF(NL) , T(MMAX.IV)

COMMON /BLKPH/ PH(MAXNl)

REAL MW

TEMP = T(K,1)

91
c

c

c

c

cc

c

cc

c

c

cc

MW = T(K,2)

TMPINF = ALF(l)

A = ALF(2)

B = ALF(3)

ALPHA = ALF(4)

PH(1) = 1.0 / ( TEMP - TMPINF ) PH(2) = ALOG10 ( A*MW + B*MW**ALPHA )

ETA = C(1)*PH(1) + PH(2)

RETURN

END

SUBROUTINE DERIV ( K, T, MMAX. IV, ALF . NL. NP1 )

PARAMETER ( MAXN1 = 11 , MAXKG = 8 )

DIMENSION ALF(NL),T(MMAX,IV)

COMMON /BLKDPH/ DPH(MAXN1 , MAXKG)

REAL MW, L0G1OE

TEMP = T(K.l)

MW = T(K,2)

TMPINF = ALF(l)

A = ALF(2)

B = ALF(3)

ALPHA = ALF (4)

92
L0G10E = ALOGIO ( EXP(l.O) ) C

DO 10 J=1,NL

DO 10 1=1, NP1

DPH(I.J) = 0.0

10 CONTINUE

C

C

C DPH(I.J) SHOULD BE SET TO THE PARTIAL DERIVATIVE

C OF PH(I) WITH RESPECT TO ALF(J)

C

DPH(l.l) = 1.0 / ( TEMP - TMPINF )**2

C

DENOM = A*MW + B*MW**ALPHA

C

DPH(2,2) = LOG10E*MW / DENOM

DPH(2.3) = LOG 10E*MW**ALPHA / DENOM

DPH(2.4) = LOG10E*B*ALOG(MW)*MW**ALPHA / DENOM

C

RETURN

END

C

C

Cc

c

c

c

SUBROUTINE WT ( M, WEIGHT, Y )

DIMENSION Y(M) .WEIGHT(M)

DO 10 1=1,

WEIGHT(I) = 1.0/Y(I)**2

10 CONTINUE

RETURN

END

93
c

c

cc

c

c

c

c

c

c

c

c

SUBROUTINE TRANS ( T, BT, Y, BY, IV, M, MMAX )

DIMENSION T(MMAX.IV), BT(MMAX.IV), Y(M), BY(M)

DO 20 1 = 1,

BY(I) = Y(I)

DO 10 J=1,IV

BT(I.J) = T(I.J)

10 CONTINUE

Y(I) = AL0G10 ( Y(I) )

20 CONTINUE

RETURN

END

SUBROUTINE BAKTRN ( ETA, BETA, M )

DIMENSION ETA(M). BETA(M)

DO 10 1=1.

BETA(I) = 10.0**ETA(I)

10 CONTINUE

C

c

c

c

94
RETURN

END

These subroutines were written into a file named VSCALC which was compiled and loaded with the INVAR programs and then executed in a terminal

session to be described in the following section.

10.3 The Terminal Session

The terminal session output that follows was obtained by executing INVAR2

on the NBS CYBER 855, running from a Qume QVT-201 terminal in screen mode. Since the terminal was not a plotting terminal, the DISSPLA plots were done after execution on the QMS laser printer with address B225.

Before reading this example, the following points should be noted:

1. Since the terminal was running in screen mode, the INVOKE,INVAR

command generated a new screen shown between the horizontal lines at the beginning of the session. This was done by the interactive NOS

Procedure INVARX in order to quiz the user for the location of the

user-supplied subroutines, the address of the QMS printer for the plots

and hardcopy output, and the name of the data file. Non-NBS users need not concern themselves with these responses.

2. The terminal session given here is not the first run of the problem.

In a previous run the logarithmic transformation of the viscosity data was performed and the resulting transformed data was written into a

file called VSTRDAT. This file is used as input for the current run so

it is not necessary to perform the transformations again.

3. The file VSTRDAT also contained the final parameter estimates from

the previous run. One of these was altered and the resulting set was

used as initial estimates for the current run.

4. After each line printer graph, the program stops with a question mark

prompt (?) generated by the TK4014 program used at NBS to generate

DISSPLA plots for the laser printer. We replied to this prompt with a "C followed by a carriage return since a carriage return alone would

have caused an end-of-file error stop.

The terminal session follows:

/INVOKE .INVAR

95
INTERACTIVE NONLINEAR LEAST SQUARES PROGRAM

CALCF FORTRAN file of user-supplied subroutines: VSCALC_

QMSADDR QMS Address for output (Default = NONE) : B22B DATAF Data file (Default = NONE) : VSTRDAT

INV102 INVAR1 or INVAR2 (Default = INVAR2)

:

Specify values and press NEXT when ready INVARX

PROC , INVARX , CALCF [ , QMSADDR] [ , DATAF] [ . INV102]

.

Executes one of two versions, INVAR1 or INVAR2, of INVAR, an interactive, variable separable, nonlinear least squares program.

Local file CALCF contains FORTRAN coded, user-supplied subroutines.

Hardcopy output and plots can be sent to QMS laser printer with

address QMSADDR.

Permanent file DATAF contains data to be fit (not needed if data

already in a local file, or is to be typed in during execution).

F3 HELP F5 HELP F6 QUIT

ENTER TITLE (MAXIMUM OF 40 CHARACTERS)

:

? VISCOSITY STUDY

DO YOU WANT DISSPLA PLOTS? (TYPE "Y" FOR YES)

? Y

DO YOU WANT THESE PLOTS SHOWN INTERACTIVELY ON

A PLOTTING SCREEN (TYPE "1"); OR SAVED IN A FILE TO BE PRINTED OUT LATER (TYPE "2")?

? 2

ENTER LABEL FOR INDEPENDENT VARIABLE 1 (MAX. OF 47 CHARACTERS):

? TEMPERATURE

ENTER LABEL FOR DEPENDENT VARIABLE (MAX. OF 47 CHARACTERS):

? LOG VISCOSITY

DO YOU WANT LINE PRINTER OUTPUT?

? Y

ARE YOUR DATA, WEIGHTS, INITIAL VALUES, ETC. ALREADY IN A VARPRO-DEFINED FILE ?

96
? Y

TYPE NAME OF FILE (MAX. OF 26 CHARACTERS)

? VSTRDAT

WOULD YOU LIKE YOUR INITIAL VALUES FILE PRINTED?

? Y

NO. OF DATA POINTS: 159

NO. OF LINEAR PARAMETERS: 1

NO . OF NONLINEAR PARAMETERS

NO. OF INDEPENDENT VARIABLES:

NDEX WEIGHT DEP VAR IND VAR 1 IND VAR 2

1 1 OOOOOOE+00 1 039110E+01 1 034000E+02 1 939992E+04

2 1 OOOOOOE+00 1 042550E+01 1 034000E+02 1 939992E+04

3 1 OOOOOOE+00 1 157900E+01 1 044000E+02 1 059986E+05

4 1 OOOOOOE+00 1 111420E+01 1 055000E+02 1 059986E+05

5 1 OOOOOOE+00 9 325700E+00 1 092000E+02 4 400480E+04

6 1 OOOOOOE+00 9 458500E+00 1 092000E+02 4 400480E+04

7 1 OOOOOOE+00 8 881800E+00 1 101000E+02 1 939992E+04

8 1 OOOOOOE+00 9 621300E+00 1 137000E+02 1 060010E+05

9 1 OOOOOOE+00 8 134600E+00 1 140000E+02 1 939992E+04

10 1 OOOOOOE+00 1 133620E+01 1 . 144000E+02 1 815098E+05

11 1 .OOOOOOE+00 7 444100E+00 1 . 185000E+02 1 939992E+04

12 1 . OOOOOOE+OO 7 449300E+00 1 . 185000E+02 1 939992E+04

13 1 .OOOOOOE+00 6 .785000E+00 1 . 200000E+02 1 .250259E+04

14 1 .OOOOOOE+OO 7 . 829000E+00 1 .201000E+02 5 .046613E+04

15 1 . OOOOOOE+00 7 .656400E+00 1 . 206000E+02 4 .400480E+04

145 1 OOOOOOE+OO 1 792000E+00 1 959000E+02 1 250259E+04

146 1 OOOOOOE+00 5 825600E+00 1 992000E+02 3 900318E+05

147 1 OOOOOOE+00 5 884100E+00 1 993000E+02 3 900318E+05

148 1 OOOOOOE+00 4 238000E+00 2 002000E+02 1 169499E+05

149 1 OOOOOOE+OO 1 562000E+00 2 005000E+02 1 109175E+04

97
150 1 000000E+00 2 227000E+00

161 1 000000E+00 4 144600E+00

162 1 000000E+00 2 392000E+00

163 1 . 000000E+00 3 307000E+00

164 1 . 000000E+00 3 849000E+00

166 1 000000E+00 3 828000E+00

166 1 . 000000E+00 1 988000E+00

167 1 . 000000E+00 3 140000E+00

168 1 . 00O000E+00 2 976000E+00

169 1 . 000000E+00 1 695000E+00

2 . 009000E+02

2.016000E+02

2 . 024000E+02

2 . 040000E+02

2.088000E+02

2.091000E+02

2 . 096000E+02

2.114000E+02

2 . 176000E+02

2 . 222000E+02

2 . 697739E+04

1.069986E+06

3.780071E+04

7.943282E+04

1 . 169499E+06

1.069986E+06

2 . 697739E+04

7.943282E+04

7.943282E+04

2 . 697739E+04

NONLINEAR PARAMETERS:

1 4.784193E+01

2 9.380365E-08

3 3.626116E-19

4 3.528796E+00

NUMERICAL DERIVATIVE INCREMENTS:

1 . OOOOOOE+OO

2 . OOOOOOE+OO

3 . OOOOOOE+OO

4 . OOOOOOE+OO

TERMS WHICH DEPEND ON NONLIN PARM 1: 10

TERMS WHICH DEPEND ON NONLIN PARM 2: 2

TERMS WHICH DEPEND ON NONLIN PARM 3: 2

TERMS WHICH DEPEND ON NONLIN PARM 4: 2

DOES DATA FOLLOW AN EVENLY SPACED TIME SERIES?

? N

NOTE: THE COMMAND "SC" WILL SAVE YOUR DATA FILE AT ANY TIME

NEXT STEP? (TYPE "H" FOR HELP)

? CALF

WHICH NONLINEAR PARAMETER IS TO BE CHANGED?

98
(TYPE ITS INDEX NO. AS IT APPEARS IN THE LIST)

? 3

ENTER NEW VALUE FOR PARAMETER 3:

? 1.0E-20

DO YOU WANT TO CHANGE ANY OTHERS?

? N

NEXT STEP? (TYPE "H" FOR HELP)

? RV

DO YOU WISH TO CHANGE ANY DEFAULTS IN THE LEAST SQUARES PROGRAM?

(BO: MAX. NO. OF ITERATIONS)

(l.OE-B: CONVERG. TOLERANCE)

(1: OUTPUT PRINTING INDEX)

? N

1 VISCOSITY STUDY

WEIGHTED NORM OF RESIDUAL = 7 .8237B56E+00

NU = .lOOOOOOE+01

ITERATION 1 NONLINEAR PARAMETERS

4.4009180E+01 3.7498487E-08

1.76B31B8E-20 3.6340762E+00

1 WEIGHTED NORM OF RESIDUAL = 2 . 8B93973E+00

NU = . BOOOOOOE+OO

NORM ( DELTA -ALF) / NORM(ALF) = 8.683E-02

ITERATION 2 NONLINEAR PARAMETERS

4.2417711E+01 3.266B216E-08

2.166B49BE-20 3 .6639208E+00

1 WEIGHTED NORM OF RESIDUAL = 2 .0880991E+00

NU = .2600000E+00

NORM(DELTA-ALF) / NORM(ALF) = 3.739E-02

99
ITERATION 3 NONLINEAR PARAMETERS

4.3127340E+01 3 . 8470079E-08

2.5B42614E-20 3.6717600E+00

1 WEIGHTED NORM OF RESIDUAL 2.0585942E+00

NU = .12BOOOOE+00

NORM(DELTA-ALF) / NORM(ALF) = 1.640E-02

ITERATION 4 NONLINEAR PARAMETERS

4.B077B62E+01 6 .3682116E-08

3.7296374E-20 3.6666329E+00

1 WEIGHTED NORM OF RESIDUAL = 2 . 0383975E+00

NU = .62BOOOOE-01

NORM (DELTA-ALF) / NORM(ALF) = 4.312E-02

ITERATION 17 NONLINEAR PARAMETERS

4.7841930E+01 9.3803727E-08

3.6260423E-19 3. B287977E+00

1 WEIGHTED NORM OF RESIDUAL = 1 . 9930702E+00

NU = 1647949E-02

NORM(DELTA-ALF) / NORM(ALF) = 1.424E-06

ITERATION 18 NONLINEAR PARAMETERS

4.7841931E+01 9.38036B4E-08

3.62B1113E-19 3.6287962E+00

1 WEIGHTED NORM OF RESIDUAL = 1 . 9930702E+00

NU = .8239746E-03

NORM(DELTA-ALF) / NORM(ALF) = 3.974E-08

100
LINEAR PARAMETERS

.7046766E+03

NONLINEAR PARAMETERS

.4784193E+02 .938036BE-07

WEIGHTED NORM OF RESIDUAL =

WEIGHTED ESTIMATED VARIANCE

.36251 HE- 18 .3628796E+01

.1993070E+01

.2B79434E-01

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED

BY TYPING "Y": TABLE OF OBSERVED * PREDICTED VALUES, WTD RESIDUAL:

? Y

COVARIANCE MATRIX:

? N

( ? Y

1

CORRELATION MATRIX:

NDEX IND VAR 1 OBS VAL PRED VAL WTD RES ID WEIGHTS

1 1 03400E+02 1 03911E+01 1 00472E+01 3 43947E-01 1 OOOOOE+00

2 1 03400E+02 1 04255E+01 1 00472E+01 3 78347E-01 1 OOOOOE+00

3 1 04400E+02 1 15790E+01 1 17734E+01 -1 94402E-01 1 OOOOOE+00

4 1 05500E+02 1 11142E+01 1 15357E+01 -4 21502E-01 1 OOOOOE+00

5 1 09200E+02 9 32570E+00 9 59673E+00 -2 71031E-01 1 OOOOOE+OO

6 1 09200E+02 9 45850E+00 9 59673E+00 -1 38231E-01 1 OOOOOE+00

7 1 10100E+02 8 88180E+00 8 68219E+00 1 99613E-01 1 OOOOOE+00

8 1 13700E+02 9 62130E+00 1 00140E+01 -3 92717E-01 1 OOOOOE+00

9 1 14000E+02 8 13460E+00 8 01496E+00 1 19644E-01 1 .OOOOOE+OO

10 1 14400E+02 1 13362E+01 1 07099E+01 6 26285E-01 1 .OOOOOE+00

11 1 18500E+02 7 44410E+00 7 33660E+00 1 07500E-01 1 . OOOOOE+OO

12 1 . 18500E+02 7 . 44930E+00 7 33660E+00 1 12700E-01 1 . OOOOOE+00

13 1 . 20000E+02 6 .78500E+00 6 .87183E+00 -8 68340E-02 1 . OOOOOE+OO

14 1 . 20100E+02 7 . 82900E+00 8 .03168E+00 -2 02682E-01 1 . OOOOOE+00

15 1 . 20600E+02 7 .65640E+00 7 .79727E+00 -1 .40873E-01 1 . OOOOOE+OO

101
145 1. 95900E+02 1. 79200E+00 1. 86556E+00 -7. 35603E-02 1 . OOOOOE+OO

146 1. 99200E+02 5 82660E+00 5 94568E+00 -1. 20081E-01 1 . OOOOOE+OO

147 1 99300E+02 5 88410E+00 5 94261E+00 -5 85071E-02 1 . OOOOOE+OO

148 2 00200E+02 4 23800E+00 4 08523E+00 1 52767E-01 1 . OOOOOE+OO

149 2 00500E+02 1 56200E+00 1 66081E+00 -9 88084E-02 1 . OOOOOE+OO

1B0 2 00900E+02 2 22700E+00 2 . 21664E+00 1 03581E-02 1 . OOOOOE+OO

151 2 O15O0E+O2 4 . 14460E+00 3 . 90006E+00 2 44540E-01 1 . OOOOOE+OO

152 2 02400E+02 2 39200E+00 2 .49895E+00 -1 06954E-01 1 . OOOOOE+OO

153 2 04000E+02 3 30700E+00 3 . 40643E+00 -9 94273E-02 1 . OOOOOE+OO

154 2 08800E+02 3 . 84900E+00 3 .83811E+00 1 .08884E-02 1 . OOOOOE+OO

155 2 09100E+02 3 . 82800E+00 3 68392E+00 1 .44076E-01 1 . OOOOOE+OO

156 2 09600E+02 1 .98800E+00 1 .96902E+00 1 .89788E-02 1 . OOOOOE+OO

157 2 . 11400E+02 3 . 14000E+00 3 . 20226E+00 -6 22605E-02 1 . OOOOOE+OO

158 2 . 17600E+02 2 .97500E+00 3 .04491E+00 -6 99061E-02 1 .OOOOOE+OO

159 2 22200E+02 1 .69500E+00 1 .66421E+00 4 .07915E-02 1 . OOOOOE+OO

CORRELATION MATRIX

ROW COLUMN COLUMN COLUMN COLUMN

1 l.OOOOOOE+OO -9 877161E-01 -9 708057E-01

2 -9.877161E-01 1 OOOOOOE+OO 9 367497E-01

3 -9.708057E-01 9 367497E-01 1 OOOOOOE+OO

4 -4.157271E-01 3 886600E-01 3 334069E-01

5

)W

-4.404127E-04

COLUMN 5

1 502759E-02 8 473789E-02

1 -4.404127E-04

2 1.502759E-02

3 8.473789E-02

4 -9.063612E-01

5 l.OOOOOOE+OO

-4.157271E-01

3.886600E-01

3.334069E-01

1 . OOOOOOE+OO

-9.063612E-01

PLOT ROUTINE CURRENTLY SET NOT TO DO ANY PLOTS.

DO YOU WANT TO START PLOTTING?

? Y

WHICH INDEPENDENT VARIABLE SHOULD BE USED?

(TYPE ITS NUMBER)

? 1

102
IF YOU WANT TO PLOT RESIDUALS, TYPE "1"

IF YOU WANT TO PLOT OBSERVED * PREDICTED VALUES, TYPE "2"

IF YOU WANT TO PLOT BOTH, TYPE "3":

? 3

RESIDUALS

IND VAR 1

1 Y DIV= .20956E-01 AND 1 X DIV= 2.3760

_j ! ! 1 j_

.6263 - *

1-

.4167

.2072

-.0024

** *

*

* * * * *

* * *

* * * 2

2 * * *3 * * * ****** * *2* ** * 2

* * * * * * *** 2 * * * **

* * * 2 * * * ** 2*2* * * 2 * 2 * * * *2 22 * * * *** * *

* ******

* * * * 2* *2

* 2*2****

2 * * *

103
I * -.2119

? C

1

.4215 - *

I _j ! ! j ! j_

103.4000 127.1600 150.9200 174.6800 198.4400 222.2000

OBSERVED (*) ft PREDICTED ( + ) VALUES

IND VAR 1

1 Y DIV= .20423 AND 1 X DIV= 2.3760

_j 1 1 1 x

_ 1- 11.7734

9.7311

7.6888

3 * **

2

+

2 +

23

+ 2*2

*

* +2

2 2*+ 2 2

+ + 22

2 * + 2

104
6 . 6466

3 . 6043

1.5620 -

I

4 ** 2 2

+ 2224 + *

* 22 2 * 3 4

2 +32

2 * 2 + * *22

42 2 2 * + +32

2 * + 22 *2+ * 2 2 * + + 2 2 *

* 22 2 2 2 2 + 2* 2 2 2

+ 2 2 +

2 2 22 + + 3

+ * * +

* + 3 2 2 2

* 2 : 4 * 2 2 3 + 2 2 * + 2 * + + 2

* *+ : * 2

2

_j j j j ! j103.4000 127.1600 150.9200 174.6800 198.4400 222.2000

? C

DO YOU WANT TO CONTINUE PLOTTING USING A DIFFERENT

INDEPENDENT VARIABLE?

? Y

WHICH INDEPENDENT VARIABLE SHOULD BE USED?

(TYPE ITS NUMBER)

? 2

TYPE IN ITS LABEL (MAX. OF 47 CHARACTERS)

? MOLECULAR WEIGHT

IF YOU WANT TO PLOT RESIDUALS, TYPE "1"

IF YOU WANT TO PLOT OBSERVED * PREDICTED VALUES, TYPE "2"

IF YOU WANT TO PLOT BOTH, TYPE "3":

? 3

1 RESIDUALS

IND VAR 2

1 Y DIV= .20956E-01 AND 1 X DIV= 7800.6

105
_ x 1 J j J j. .6263 - * I

I

I

I

I

II * I

I .4167 - 2

I

I * I * I

II 2

I13 *

I * .2072 - * I

I 5

I * * *

I &

I X

I 7

I 9

16 * 18 2 * .0024 - 8

I X

I 6 1X22 *

I 7 16 3 *

I 8 * *

I 6

II 3

.2119 - *

I * I 2

I 3

I 2

II

106
I * I

I * .4216 - *

? C

_! j J J ! !_

.OOOE+00 .780E+0B .156E+06 . 234E+06 .312E+06 .390E+06

OBSERVED (*) t PREDICTED (+) VALUES

IND VAR 2

1 Y DIV= .20423 AND 1 X DIV= 7800.6

_j J ! j ! !_

11.7734 - +

I 2

I * I * I +

I 2

I

I 2

I 2

9.7311

I +

- +

I 6

I 4

I 2

I 4

I 4

I 2

I 7

I 6

I 6

7.6888 - 6

I 5

I 7

I 7

I X

I X

I X

I X

I X

I X 6.6466 - 8

I X

I X

I 4

I X

107
I 7

I X

I 7 * *

I 5 +

I 8 2 * + + 3.6043 - 2 + * -

I 9 2

I 3 2

I 9 2

I 3 +

I 4 2

I 2 *

I 2 2

I 2 + 2

I * 3 1.6620 - 3

-

I -I

.001 —I—

.780E . 156E+06

---I .234E+06

1 .312E+061-

DEt 00 +05 . 390E+06

? c

DO YOU WANT TO CONTINUE PLOTTING USING A DIFFERENT

INDEPENDENT VARIABLE?

? N

NEXT STEP? (TYPE "H" FOR HELP)

? RS

INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED: PARAMETERS, THEIR STD DEV ft T-RATIO:

? Y

RES SUM OF SqS, RES MEAN SQ , RES STD DEV ft COEFF OF DETERMINATION:

? Y

TABLE OF OBS. PRED VAL, RESID ft STD RESID:

? Y

PLOT OF STD RES VS IND VAR:

? Y

PLOT OF STD RES VS PRED VAL:

? Y

NORMAL PROBABILITY PLOT:

? Y

108
PARAMETER S.D. OF PARM. T-RATIO

7 . 046766E+02

4.784193E+01

9 . 380365E-08

3.625111E-19

3.628796E+00

3.001868E+01

1.903968E+00

3.B14149E-08

3.168989E-19

6. 666 11 BE -02

347460E+01

612761E+01

669313E+00

147BB4E+00

293632E+01

WTD RESIDUAL SUN OF SQUARES: 3.972329E+00

WTD RESIDUAL MEAN SQUARE: 2.679434E-02

WTD RESIDUAL STANDARD ERROR: 1.606062E-01

COEFFICIENT OF DETERMINATION (R-SQUARE) 9.947714E-01

INDEX IND VAR 1 DBS VAL PRED VAL RES ID STD RES

1 1.03400E+02

2 1 . 03400E+02

3 1.04400E+02

4 1 . 06600E+02

6 1 . 09200E+02

6 1.09200E+02

7 1 . 10100E+02

8 1 . 13700E+02

9 1 . 14000E+02

10 1 . 14400E+02

11 1.18500E+02

12 1.18600E+02

13 1 . 20000E+02

14 1 . 20100E+02

15 1.20600E+02

1.03911E+01

1.04265E+01

1.15790E+01

1.11142E+01

9.32570E+00

9 . 458B0E+00

8.88180E+00

9.62130E+00

8.13460E+00

1 . 13362E+01

7.44410E+00

7 . 44930E+00

6.78600E+00

7 . 82900E+00

7.65640E+00

1.00472E+01

1.00472E+01

1.17734E+01

1.16367E+01

9 . B9673E+00

9.B9673E+00

8.68219E+00

1.00140E+01

8.01496E+00

1.07099E+01

7 . 33660E+00

7 . 33660E+00

6.87183E+00

8.03168E+00

7.79727E+00

3.43947E-01

3.78347E-01

-1.94402E-01

-4.21B02E-01

-2.71031E-01

-1.38231E-01

1.99613E-01

-3.92717E-01

1.19644E-01

6.26286E-01

1.07BOOE-01

1.12700E-01

-8.68340E-02

-2.02682E-01

-1.40873E-01

7.73207E-02

7.99B18E-01

-4.14347E-01

-9.08277E-01

-6.07B69E-01

-3.09873E-01

4.B2273E-01

-9.26B3BE-01

2.83BB1E-01

1.48918E+00

2.69062E-01

2.82066E-01

-2.21284E-01

-B.16977E-01

-3.616B6E-01

146 1 . 96900E+02

146 1 . 99200E+02

147 1 . 99300E+02

148 2 . 00200E+02

1.79200E+00 1.86BB6E+00 -7.36603E-02 -7.096B2E-01

6.82560E+00 B .94668E+00 -1 . 20081E-01 -1.2B8BBE+00

6.88410E+00 6.94261E+00 -B.8B071E-02 -6.15B07E-01

4.23800E+00 4.08B23E+00 1.62767E-01 1.76382E+00

109
149 2 00500E+02 1 66200E+00 1 66081E+00 -9 88084E-02 -1 13245E+00

160 2 00900E+02 2 22700E+00 2 21664E+00 1 03581E-02 1 23628E-01

151 2 01600E+02 4 14460E+00 3 90006E+00 2 44540E-01 2 99568E+00

152 2 02400E+02 2 39200E+00 2 49895E+00 -1 06954E-01 -1 36252E+00

163 2 04000E+02 3 30700E+00 3 40643E+00 -9 94273E-02 -1 38257E+00

154 2 08800E+02 3 84900E+00 3 83811E+00 1 08884E-02 2 29712E-01

155 2 09100E+02 3 82800E+00 3 68392E+00 1 44076E-01 3 14602E+00

156 2 09600E+02 1 98800E+00 1 96902E+00 1 89788E-02 4 47527E-01

167 2 11400E+02 3 14000E+00 3 20226E+00 -6 22606E-02 -2 03067E+00

168 2 17600E+02 2 .97600E+00 3 04491E+00 -6 99061E-02 -1 50806E+00

159 2 22200E+02 1 69500E+00 1 .66421E+00 4 07915E-02 6 .21766E-01

ENTER

? 1 ANY INTEGER ' m CONTINUE

STANDARDIZED RESIDUAL VS INDEPENDENT VARIABLE 1

1 Y DIV= .10983

-j J AND 1 X DIV= 2.3760

.j j j_ —I- 3.1460

2 . 0477

.9495

* **

* *

* *

* * * * *

* * * *

* * * *2

* 2 ****2** * 2 *2* * *** * * * * * *** 2 * *

110
I ++* * *2 *2* 3 * *

- . 1488 - 2 2***** *

* * * *2 2*

* 2 * * *** * * * * * * *

* * * * 2 * *** *

* * * **

* * * * * * 2 2

* **

* * 1.2471 - * * 2

-2.3464

? C

1

_ x J ! j J !_

103.4000 127.1600 150.9200 174.6800 198.4400 222.2000

STANDARDIZED RESIDUAL VS INDEPENDENT VARIABLE 2

1 Y DIV= .10983 AND 1 X DIV= 7800.6

_! J ! j !-._

3.1460 - *

I * I

I

I

I

I

I

II

1-

2 . 0477

I *

I 2 12*

I 2

I * I 3

111
I

II .9495 -

I 2

I 2

I 6

I 4

I 6

I X

I X *

I 9 *

I X -.1488 - X

I 9

I X

I 5

I X

I 3

I 3

I 7

I * 2

I * -1.2471 -1*2

I

II

I

I

II

I -2.3464 - *

_j j ! 1 j !_

.OOOE+00 .780E+05 .156E+06 . 234E+06 .312E+06 .390E+06

STANDARDIZED RESIDUAL VS. PREDICTED VALUES

1 Y DIV= .10983 AND 1 X DIV= .20238

_! ! j J J !- 3.1460 - *

I *

I

I

I

112
2 . 0477

**

.9495

-.1488

-1.2471

* * **

* * * ** *

* * * *

* ******

2 2   2 2 ** 

* 2 * * 2 *** * * ********* *

* * 2 22*** **2*

* *** ** 2 * *

* 2 * *2 **

* 2 * * * * *

* * * * * * 2 * * * *2 **

* * * * *

* * * 2 2 * *

* * *

* *

*

* * *

*

-2.3454

-I 1- I— -I 1-

113
? c

1

1.6B42 3.6780 6.7019 7.7257 9.7496 11.7734

NORMAL PROBABILITY PLOT

1 Y DIV= .10983 AND 1 X DIV= .10564

_! ! J J !-. 3 . 1460

2 . 0477

.9495

1488

1.2471

1-

* -

2*

* 2*

2 ** 33

4 24* 563

57 27*

574

262

45

262

23

353

4*

3 33*

*2

* *

114
***

-2.3454 - *

I _! ! ! ! j j_

-2.6262 -1.5698 -.5134 .5430 1.5994 2.6558

? C

NEXT STEP? (TYPE "H" FOR HELP)

? SO

END OF DISSPLA 9.0 -- 17196 VECTORS IN 8 PLOTS. RUN ON 5/8/87 USING SERIAL NUMBER AT NBS/CS**2 - GAITHERSBERG,

MD. PROPRIETARY SOFTWARE PRODUCT OF ISSCO, SAN DIEGO, CA. 4553 VIRTUAL STORAGE REFERENCES; 4 READS; WRITES.

10.4 The Graphical Output

The 8 DISSPLA plots generated by the run are shown in Figures 10.1-10.8.

Since the fitted function defines a 2-dimensional surface, the plots of data

and fit vs. the two variables T and M give projections of this surface. For

the sake of simplicity, the program does not try to plot a family of curves

to represent this surface, but rather evaluates the fitting function only at the observation points and plots both the observations and fits as pairs of

discrete points. The various plots of the residuals are not the kind of scatter diagrams that one would hope to obtain from a really good fit, and the

normal probability plot (Figure 10.8) indicates significant deviations from

the diagonal distribution that one would expect if the residuals were truly

normally distributed. It might be possible to remedy some or all of these

problems by a proper weighting of the observations. If not, then one would

be forced to conclude that the model is not quite right.

115
RESIOURLS

<6 "

K "

IT) _)

2

^ -

,<v

° 3

oOq

Po.

° a °

8 B o o o oo S

O D o **Tftoqf«o° » b-t»

o J

5 *o a« oo .IobdO o o

o o 1° c?

N

#*

N

#' £>

TEMPERATURE

,f>'

T

,•?'

Figure 10.1: Residuals (data- model) vs. temperature

116
VISCOSITY STUDY

OBSERVED - , PREOICTEO -

<y

N<>'

tno

CDO Q

_J K'

o -

o

6°

9

9 o _~ •

•• ••* • 9 « O

• • * 69

t

.«$>•T

s#'

T

s# <* ^'

TEMPERATURE

Figure 10.2: Data and fit vs. temperature

117
RESIDUALS

: MOLECULAR WEIGHT "10

Figure 10.3: Residuals (data - model) vs. molecular weight

118
VISCOSITY STUDY

OBSERVED - , PREDICTED -

: MOLECULAR WEIGHT

Figure 10.4: Data and fit vs. molecular weight

119
STRND. RESID. VS. IND. VRR.

CE 7

IT) ^• u

O N

Ld

ES *

(E

I—

tn <>> - ,V

tfl

° o - o o o

o °0

1 Z n

•a B

•"

o

o

#*

^

o

s#- <«>'

6 o D o i

«P„

O ^ O

#' s*' #'

TEMPERATURE

#'

Figure 10.5: Standardized residuals vs. temperature

120
STAND. RESID. VS. IND. VRR.

: MOLECULAR WEIGHT

Figure 10.6: Standardized residuals vs. molecular weight

121
STRND. RESIDURLS VS. PRED. VALUES

3^

2 o 1 en V

LJa v

LJ

oc

<nIN

'

en ^ - ,v

o<ft

IP^>°A^£°^>

o°(P*o

- »

o

a o a " « A•a u

ao

T T

K' <o- *>' LOG (VISCOSITY)

s*'

o a

<*

Figure 10.7: Standardized residuals vs. viscosity predicted by fit

122
NORMRL PROBRBILITY PLOT

A* A* <V v V

EXPECTED STflNDRRDIZED RESIDURLS

Figure 10.8: Standardized residuals versus the expected standardized residuals when it is assumed that the measurement errors are independently

normally distributed

123
Appendix A

THE SUBROUTINES IN

THE INVAR CODE

Figure A.l depicts the INVAR subroutine network. The arrows represent

the directions of the subroutine CALLs. The five user supplied subroutines

are written in boldface. They are described in detail in Chapter 3. The

other subroutines are: TYPE The subroutine which nominates the DISSPLA plotting devices.

PTEKAL nominates Tektronics plotting terminals. TK4014 and

HWSPEC nominate the QMS laser printers.

PR The subroutine used to print out the various values requested by the

user. CH The subroutine used to change various quantities at the user's request.

ARRAYD A subroutine use to print out matrices.

VARPRO The main driver for the VARPRO subroutine package. The

. VARPRO subroutines ORFACl, ORFAC2, DPA, BACSUB, POSTPR,

COV, and VARERR are unchanged, but ADA was altered and

NUMBER was added by Rail and Funderlic [12]. PLOT1, PLOT2 Subroutines used to generate the DISSPLA plots.

STAT The subroutine which generates the various statistical diagnostics

obtainable with the "RS" command.

ORPLOT A sorting subroutine used by STAT.

VPLT, VPLT2, VMAINH, VPLTER, VPMAIN Line printer plot routines.

124
MAIN

r»PTEKAL

*TYPE-4-»TK4014

TRANS L»HWSPEC

PR

WT

TCH

ARRAYD

VARPROBAKTRN

- PLOT 2-

-•ORFAC2

-POSTPR1 -ORFAC1

i

-»VARERR

t — DPA

1 + BACSUB

COV

DEMV«-

-» ADA

CALC«-

i > disspla*-

-*PLOT

T —• STAT

-»VPLT

Jl.

P

ORPLOT

VPLT2-

VMAINH

»VPLTER

-VPMAIN

-»NUMDER

TYPEND-

-•SETUP2-

SPLT-

-SETUP1 i_

-»PLTMN<

-PLTERR

SPLTH

it

DONEPL

Figure A.l: The INVAR subroutine hierarchy

125
SETUP1, SETUP2, SPLT, SPLTH, PLTMN, PLTERR Line printer

plot routines.

TYPEND An entry in Subroutine TYPE which is called after the user enters any of the stop commands. It terminates the DISSPLA plotting

by calling the DISSPLA subroutine DONEPL.

126
Appendix B

CHANGING THE INVAR

PROGRAMS

B.l Changing the PROGRAM Statement

The PROGRAM statement in the MAIN program has the form

PROGRAM MAIN (INPUT ,OUTPUT=/1000 .LPOUT , INVRPL=/1000,

C

1 TAPE5=INPUT,TAPE6=0UTPUT,TAPE10=LP0UT)

which is non-standard FORTRAN specific for Cyber computers. User who

want to run INVAR on some other computer will want to replace or remove

this statement. TAPE5 and TAPE6 are just the standard input and output

units. INVRPL is the file for the DISSPLA plots, and TAPE10 is the unit

for the hardcopy line printer output. User who change this PROGRAM

statement must do whatever is necessary to assign, open and close these

files on their own computer.

B.2 Changing the DISSPLA Plotting Device Nornination

The DISSPLA plotting device nomination is done in Subroutine TYPE.

Currently there are two ways to get the plots: in real-time on a Tektronics

plotting terminal or after execution has been completed on a QMS laser printer.

127
To change the real-time nomination to some device other than a Tektronics terminal, change

IF (IDP) CALL PTEKAL

to

IF (IDP) CALL devnam

where "devnam" is the appropritate device name. To change the deferred, hard-copy plotter nomination to some device

other than a QMS laser printer, change

IF (IDCP) THEN

CALL TK4014 (960.1)

CALL HWSPEC ( ' INVRPL' , 'FILE'

)

END IF

to

IF (IDCP) THEN

CALL devnam

END IF

where "devnam" is the appropriate device name.

B.3 Changing the Array Dimensions

Some problems may require larger arrays than are presently available in IN- VAR. Arrays can easily be redimensioned by changing a few PARAMETER

statements. Each of the subsections in this section lists, by subroutine, the

alterations required to change one of the dimensioning parameters.

Changing the maximum number of observations

To change the maximum number of observations from 300 to m:

• In PLOT1 replace

PARAMETER (MMAX=300)

with

128
PARAMETER (MMAX= m )

• In MAIN replace

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=33,

\ MAXN1=11, MAXNKG=18. MAXKG=8 . MAXP=21)

with

PARAMETER (MMAX= m , MMAX2=2*m, MIV=6, MAXN=10, MAXNP2=33,

2 MAXN1=11, MAXNKG=18, MAXKG=8, MAXP=21)

• In ORPLOT replace

PARAMETER (MMAX=300)

with

PARAMETER (MMAX= m )

• In POSTPR replace

PARAMETER (MMAX=300, MAXNKG=18)

with

PARAMETER (MMAX= m , MAXNKG=18)

• In STAT replace

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=18, MIV=6)

with

PARAMETER (MMAX= m , MAXN1=11. MAXNKG=18, MIV=6)

129
Changing the maximum number of independent variables

To change the maximum number of independent variables per observation

from 6 to iv:

• In subroutine CH, the format to print the data points (FORMAT

statement 60) will need to be changed.

• In MAIN replace

PARAMETER (MMAX=300. MMAX2=600, MIV=6, MAXN=10, MAXNP2=33,

2 MAXN1=11. MAXNKG=18, MAXKG=8 , MAXP=21)

1650 FORMAT (7(E14.8,2X) )

with

PARAMETER (MMAX=300. MMAX2=600, MIV=iv, MAXN=10, MAXNP2=33

MAXN1=11, MAXNKG=18, MAXKG=8, MAXP=21)

1650 FORMAT ( (iv+1) (E14.8.2X) )

• In subroutine PR, the formats to print out data points and weights

(FORMAT statements 410, 450 and 480) will need to be changed.

• In STAT replace

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=18, MIV=6)

with

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=18, MIV=iv)

Changing the maximum number of linear parameters

To change the maximum number of linear parameters from 10 to n (k =

max. no. of nonlinear parameters; p = max. no. of nonzero partials):

130
• In ADA replace

PARAMETER (MAXN1=11, MAXN=10, MAXKG=8 . MAXP=21)

with

PARAMETER (MAXNl=n+l, MAXN=n, MAXKG=8, MAXP=21)

• In CALC replace

PARAMETER (MAXN1=11)

with

PARAMETER (MAXN=n+l)

• In DERIV replace

PARAMETER (MAXN1=11, MAXKG=8)

with

PARAMETER (MAXNl=n+l, MAXKG=8)

• In DPA replace

PARAMETER (MAXN1=11, MAXKG=8)

with

PARAMETER (MAXNl=n+l , MAXKG=8)

• In MAIN replace

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=33

I MAXN1=11, MAKNKG=18, MAXKG=8 , MAXP=21)

710 FORMAT (38A1.11I2)

131
1750 FORMAT (' TERMS WHICH DEPEND ... ',12,''. \11I2)

with

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=n, MAXNP2=n+p+2

,

2 MAXNl=n+l, MAXNKG=n+k, MAXKG=8. MAXP=21)

710 FORMAT (38A1 , (n+1) 12)

1750 FORMAT TERMS WHICH DEPEND ... ',12,': \(n+l)I2)

• In NUMDER replace

PARAMETER (MAXN=10. MAXN1=11, MAXKG=8, MAXP=21)

with

PARAMETER (MAXN=n, MAXNl=n+l , MAXKG=8, MAXP=21)

• In POSTPR replace

PARAMETER (MMAX=300. MAXNKG=18)

with

PARAMETER (MMAX=300, MAXNKG=n+k)

• In PR replace

132
PARAMETER (MAXN1=11, MAXKG=8 , MAXP=21)

320 FORMAT (/' '/' TERMS WHICH DEPEND ...

, ,I2, ,

: ',1112)

with

PARAMETER (MAXNl=n+l, MAXKG=8, MAXP=21)

320 FORMAT (/' '/' TERMS WHICH DEPEND ... ',12,': \(n+l)I2)

• In STAT replace

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=18, MIV=6)

with

PARAMETER (MMAX=300, MAXNl=n+l, MAXNKG=n+k, MIV=6)

Changing the maximum number of nonlinear parameters

To change the maximum number of nonlinear parameters from 8 to k (n =

max. no. of linear parameters):

• In ADA replace

PARAMETER (MAXN1=11, MAXN=10, MAXKG=8 , MAXP=21)

with

PARAMETER (MAXN1=11, MAXN=10. MAXKG=k, MAXP=21)

• In DERIV replace

PARAMETER (MAXN1=11, MAXKG=8)

133
with

PARAMETER (MAXN1=11. MAXKG=k)

• In DPA replace

PARAMETER (MAXN1=11, MAXKG=8)

with

PARAMETER (MAXN1=11, MAXKG=k)

• In MAIN replace

PARAMETER (MMAX=300, MMAX2=600. MIV=6, MAXN=10. MAXNP2=33,

I MAXN1=11, MAXNKG=18. MAXKG=8 , MAXP=21)

690 FORMAT (23A1 .8E14.8)

1720 FORMAT (.' NONLINEAR PARAMETERS: ' .8E14.8)

with

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=33,

2 MAXN1=11, MAXNKG=k+n, MAXKG=k, MAXP=21)

690 FORMAT (23A1 ,kE14 .8)

134
1720 FORMAT (' NONLINEAR PARAMETERS : \kE14.8)

• In NUMDER replace

PARAMETER (MAXN=10, MAXN1M1. MAXKG=8 , MAXP=21)

with

PARAMETER (MAXN=10, MAXN1=11. MAXKG=k. MAXP=21)

• In POSTPR replace

PARAMETER (MMAX=300, MAXNKG=18)

with

PARAMETER (MMAX=300, MAXNKG=k+n)

• In PR replace

PARAMETER (MAXN=11, MAXKG=8 , MAXP=21)

with

PARAMETER (MAXN=11, MAXKG=k, MAXP=21)

• In STAT replace

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=18, MIV=6)

with

PARAMETER (MMAX=300, MAXN1=11, MAXNKG=k+n, MIV=6)

Changing the maximum number of nonzero partial derivatives

To change the maximum number of nonzero partial derivatives from 21 to p (n = max. no. of linear parameters):

135
•

In ADA replace

PARAMETER (MAXN1=11, MAXN=10, MAXKG=8 , MAXP=21)

with

PARAMETER (MAXN1=11, MAXN=10, MAXKG=8, MAXP=p)

In CH replace

PARAMETER (MAXP=21)

with

PARAMETER (MAXP=p)

In MAIN replace

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=33,

2 MAXN1=11, MAXNKG=18, MAXKG=8, MAXP=21)

720 FORMAT (34A1 .21E14 .8)

1770 FORMAT (' NUMERICAL DERIVATIVE INCREMENTS: \21E14.8)

with

PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=p+n+2

,

2 MAXN1=11, MAXNKG=18, MAXKG=8 , MAXP=p)

720 FORMAT (34A1 ,pE14.8)

136
1770 FORMAT (' NUMERICAL DERIVATIVE INCREMENTS: \pE14.8)

• In NUMDER replace

PARAMETER (MAXN=10, MAXN1=11, MAXKG=8, MAXP=21)

with

PARAMETER (MAXN=10, MAXN1=11, MAXKG=8 , MAXP=p)

• In PR replace

PARAMETER (MAXN1=11, MAXKG=8 , MAXP=21)

with

PARAMETER (MAXN1=11, MAXKG=8 , MAXP=p)

B.4 Changing the Iteration Specification Parameters

Changing the maximum number of iterations

To change the maximum number of iterations from 250 to i:

• In MAIN replace

1020 FORMAT (/' NUMBER OF ITERATIONS MUST BE GREATER 7lH , 2 'THAN AND LESS THAN OR EQUAL TO 250. '/1H , 3 'HOW MANY ITERATIONS? ') READ (5,*) ITMAX

IF ( ITMAX. LE.O. OR. ITMAX. GT. 250) GO TO 1010

with

137
1020 F0RMAT(/' NUMBER OF ITERATIONS MUST BE GREATER '/1H , 2 'THAN AND LESS THAN OR EQUAL TO i

. '/1H , 3 'HOW MANY ITERATIONS? ')

READ (5,*) ITMAX

IF (ITMAX. LE.O. OR. ITMAX. GT. i ) GO TO 1010

Changing the default number of iterations

To change the default number of iterations from 50 to q: • In MAIN replace

DATA ITMAX/50/, EPS1/1 .OE-5/. IPRINT/1/

980 FORMAT (/' ... (50:MAX. NO. OF ITERATIONS)' ... )')

with

DATA ITMAX/q/, EPS1/1 .OE-5/ , IPRINT/1/

980 FORMAT (/' ... ( q :MAX. NO. OF ITERATIONS)' ... )')

Changing the maximum value of the Marquardt parameter

To change the maximum value allowed for the Marquardt parameter from

100 to f: • In VARPRO, for the third statement following statement number 800,

replace

IF (NU.LE.100) GO TO 810

with

IF (NU.LE. f ) GO TO 810

138
Appendix C

DESCRIPTION OF

VARIABLES USED IN

INVAR

A list of variables used in INVAR follows. Each variable name and, where

applicable, its dimension specifications are printed in boldface. Following

that is a short description of the variable and, enclosed in brackets, a list of

the subroutines in which it appears.

A(300,33) On input to VARPRO, the array contains the PH(J)'s and their

derivatives. On output, the first N+KG rows and columns contain an

approximation to the (weighted) covariance matrix at the solution (the

first N rows correspond to the linear parameters, the last the last KG

to the nonlinear ones), column N+KG+1 contains the weighted residuals, A(l,N+KG+2) contains the (Euclidean) norm of the weighted

residual and A(2,N+KG+2) contains an estimate of the (weighted)

variance of the observations. [ADA, COV, DPA, MAIN, POSTPR,

STAT, VARPRO]

ACUM Variable used to aid in computation, e.g., accumulate a sum. [BACSUB, DPA, ORFAC1, ORFAC2, POSTPR, VARPRO]

AI Absolute value of the Ith diagonal element of A. [STAT]

AIMl Floating point equivalent of 1-1 (dummy subscript), used in computing values for the X axis label (terminal and/or line printer plots).

[PLTMN]

AIROW Floating point equivalent of IROW. [PLTMN]

139
AJ Absolute value of the Jth diagonal element of A. [STAT]

ALF(8) Vector of nonlinear parameters. [ADA, CALC, CH, DERIV, DPA,

MAIN, NUMDER, POSTPR, PR, STAT, VARPRO]

ALFNRM Norm of ALF. [VARPRO]

ALINE(103) Character vector of the current graph line (terminal and/or

line printer plots). [PLTMN]

ALPHA Factor used in calculating Householder reflections. [DPA, ORFAC1, ORFAC2]

ALPHAI Axis symbol T (terminal and/or line printer plots). [PLTMN,

VPMAIN]

ALPHAX Plotting symbol "X" (terminal and/or line printer plots). Used

when more than 9 points fall in a single location. [PLTMN]

ANS User response to question asking for increment values. [MAIN]

ANUMLM NUMLAB-1. [PLTMN]

ASAVE Value of ALF (I) when computing numerical derivative. [NUMDER]

AXISCH Y axis symbol to be used on the current graph line (terminal

and/or line printer plots). Symbol is either "P or "-". [PLTMN]

Al(300,10) Triangular matrix used in solving system of equations. [BACSUB]

B(300,9) Matrix which contains the derivative portion of A. [ORFAC1,

ORFAC2]

BA(l) Singly-dimensioned array to be printed as a matrix. [ARRAYD]

BB Square root of AI or AJ. [STAT]

BETA Factor used in calculating Householder reflections. [DPA, ORFAC1,

ORFAC2]

BKETA Vector of the back-transformed values of the predicted Y. [BAKTRN, MAIN]

BLANK Plotting symbol (blank) (terminal and/or line printer plots). [PLTMN,

VMAINH, VPMAIN]

BLIDCP Common block containing IDCP. [MAIN, STAT, TYPE]

BLKBM Common block containing BM. [POSTPR, STAT]

BLKDPH Common block containing DPH. [ADA, DERIV, MAIN, NUMDER]

BLKIDP Common block containing IDP. [MAIN, PLOT1, PLOT2, STAT,

TYPE]

140
BLKINC Common block containing INC. [ADA, DPA, MAIN, PR]

BLKITP Common block containing ITP. [MAIN, PR, SETUP1, SETUP2,

STAT]

BLKNDI Common block containing NDI. [ADA, CH, MAIN, NUMDER,

PR]

BLKPH Common block containing PH. [ADA, CALC, MAIN, STAT]

BLKTS Common block containing TS. [MAIN, PLOT1, PLOT2, STAT]

BLOCKP Common block containing P. [ADA, MAIN]

BLOKIP Common block containing IP. [ARRAYD, CH, MAIN, POSTPR,

PR, SETUP1, SETUP2, STAT, VARERR, VARPRO]

BM(300,18) Matrix of predictor variables used in calculating the stan- dardized residuals. [POSTPR, STAT]

BT(300,6) Array of back-transformed independent variables. BT(I,J) contains the value of the Ith observation of the Jth independent variable.

[MAIN, TRANS]

BY(300) Vector of back-transformed values of the dependent variable.

[MAIN, TRANS]

Bl N+2, used to identify the first column of A which contains derivatives

of PH(J). [VARPRO]

C(10) Vector of linear parameters. [ADA, CALC, MAIN, STAT, VARPRO]

CC(10) Vector used to identify which PH(J) is to be used in calculating

the numerical derivative. [NUMDER]

COLUMN Character variable containing "COLUMN". [ARRAYD]

CR(18,18) Correlation matrix. [MAIN]

CTSS Corrected total sum of squares. [STAT]

D Plotting symbol "+" (terminal and/or line printer plots). [VMAINH,

VPMAIN]

DALFNM Norm of the change in ALF from the previous iteration to the

present one. [VARPRO]

DELR Change in the norm of the residual from the previous iteration to the present one. [VARPRO]

DELTA Distance between any 2 of the 500 points which plot the curve of

the predicted Y (DISSPLA plots). [MAIN]

141
DELX The range of the X axis (terminal and/or line printer plots). [PLTMN]

DELY The range of the Y axis (terminal and/or line printer plots). [PLTMN]

DK Numerical derivative increment. [NUMBER]

DPH(11,8) DPH(I,J) is the derivative of PH(I) with respect to ALF(J).

[ADA, DERIV, MAIN, NUMDER]

DUMMY Dummy argument in subroutine parameter list. [STAT]

DV Derivative value to be entered in A. [NUMDER]

EPSl Absolute tolerance for the norm of the projection of the residual onto

the range of the Jacobian of the variable projection functional [MAIN,

VARPRO]

EPS2 Relative and absolute tolerance for the change in the norm of the

residual. [VARPRO]

ERRCHK Common block containing IERROR. [PLTERR, SPLT, SPLTH,

VPLT, VPLTER, VPLT2]

ERRNO Error number (absolute value of IERR). [VARERR]

ETA(300) Vector of predicted values of the dependent variable. [ADA,

BAKTRN, CALC, MAIN, SETUP2, STAT]

ETAPL(600) Vector of observed and predicted values of the dependent

variable in the order they are to be plotted (terminal and/or line

printer plots). [MAIN, SETUP2]

EX Logical variable used in INQUIRE statement to indicate if a file name

already exists or not. [MAIN]

FIRSTC First column to which Householder reflections should be applied.

[DPA]

FIRSTR First row to which Householder reflections should be applied.

[DPA]

FNAME Name of file conataining all initial information (data points, parameters, variables, etc.). [MAIN]

GNSTEP Step size for Gauss-Newton iteration. [VARPRO]

H Variable used in arranging indexing of IOR. [ORPLOT]

HN Variable used in arranging indexing of IOR. [ORPLOT]

HYPHEN Axis symbol "-" (terminal and/or line printer plots). [PLTMN]

142
I Dummy subscript. [ADA, BACSUB, CH, COV, DPA, MAIN, NUMDER,

ORFAC1, ORFAC2, ORPLOT, PLOT1, PLOT2, PLTMN, POSTPR,

PR, STAT, VMAINH, VPMAIN]

IA(ll) Vector used in determining INC matrix. [MAIN, PR]

IANS User response to interactive questions. [CH, MAIN, PR, STAT]

IB Dummy subscript. [MAIN]

IBACK Dummy subscript. [BACSUB]

IC Number of values falling outside graph bounds (terminal and/or line

printer plots). [PLTMN]

ICODE Command code chosen by the user. [CH, MAIN, PR]

ICOL Dummy subscript. [PLTMN]

ICR Flag to indicate that the correlation matrix should be printed. [MAIN]

ICV Flag to indicate that the covariance matrix should be printed. [MAIN]

IDCP Flag to indicate that DISSPLA plots should be saved in a file (META)

for later access. [MAIN, STAT, TYPE]

IDP Flag to indicate that DISSPLA plots should be printed interactively on a plotting screen. [MAIN, PLOT1, PLOT2, STAT, TYPE]

IERR Integer error flag. [MAIN, ORFAC1, VARERR, VARPRO]. The

returned value is interpreted as follows:

>0 successful convergence.

-1 terminated due to too many iterations.

-2 terminated for ill-conditioning (Marquardt parameter too large).

-4 input error in parameter M, N, KG, NPP2 or MMAX.

-5 INC matrix improperly specified.

-6 weights not nonnegative.

-7 "constant" column was computed more than once. -8 catastrophic failure-a column of matrix A has become zero. IERROR Flag to indicate type of error in graph (terminal and/or line

printer plots). [PLTERR, SPLT, SPLTH, VPLT, VPLTER, VPLT2]

II Dummy subscript. [ADA, DPA, MAIN, STAT]

IK Subscript used to identify placement of derivative in matrix A. [ADA]

IK1 Dummy subscript. [PLTMN]

143
IL Dummy subscript. [CH, PR]

IMEAN Flag to indicate whether zero line is to be plotted (terminal and/or

line printer plots). [VMAINH, VPLT, VPLT2, VPMAIN]

IN Subscript to identify which numerical derivative increment should be

used. [ADA, NUMDER]

INC(8,11) Incidence matrix for partial derivatives. INC(K,J)=1 if nonlinear parameter ALF(K) appears in PH(J). Otherwise, INC(K,J)=0.

[ADA, DPA, MAIN, PR]

INCKJ INC(K,J). [DPA]

IND Character variable containing "IND". [MAIN, PR, STAT]

INV Number of independent variable used for plotting (terminal and/or

line printer plots). [MAIN, SETUP1, SETUP2]

IOR(300) Vector containing indexing of independent variable values ar- ranged in ascending order. [ORPLOT, STAT]

IP Flag to indicate line printer output. [ARRAYD, CH, MAIN, POSTPR,

PR, SETUP1, SETUP2, STAT, VARERR, VARPRO]

IPCODE Flag to indicate, during error checking, which subroutine has

been called in plotting (terminal and/or line printer plots). [PLTERR,

SPLT, SPLTH, VPLT, VPLTER, VPLT2]

IPI Flag to indicate whether to plot residuals or predicted and observed

values. [MAIN, STAT]

IPIBLK Common block containing IPI. [MAIN, STAT]

IPIV Flag to indicate whether to plot standardized residuals versus independent variables. [STAT]

IPP Flag to indicate whether to plot probability of standardized residuals.

[STAT]

IPPV Flag to indicate whether to plot standardized residuals versus pre- dicted values. [STAT]

IPR Flag to indicate whether output is sent to terminal or to line printer.

[PLTERR, PLTMN, SPLT, SPLTH, VMAINH, VPLT, VPLTER, VPLT2,

VPMAIN]

IPRINT Integer used to control output printed by VARPRO. If IPRINT>0,

the nonlinear parameters, norm of residual and Marquardt parameter

will be printed every IPRINTth iteration. If IPRINT=0, only the

final quantities will be printed. [DPA, MAIN, ORFAC1, POSTPR,

VARERR, VARPRO]

144
IPR10 Directs line printer output into temporary file 10. [MAIN, SETUPl,

SETUP2, STAT]

IPR6 Directs line printer output to be printed on terminal. [MAIN, SETUPl,

SETUP2, STAT]

IPT Flag to indicate that t-ratios should be printed. [STAT]

IP1 1+1 (dummy subscript). [BACSUB, COV]

IQ Flag to indicate Command Segment. [MAIN]

IR Flag to indicate that residual sum of squares, residual mean square,

residual standard error and coefficient of determination should be

printed. [STAT]

IROW Line of the graph currently being printed out (terminal and/or line

printer plots). [PLTMN]

ISEL Flag to indicate what (function values and/or partial derivatives)

should be computed. [ADA, DPA]

ISS Integer value of the specified symbol to be plotted (terminal and/or

line printer plots). [PLTMN]

ISTOP Element in ALINE vector which contains the last significant symbol to be printed (terminal and/or line printer plots). [PLTMN]

ISUB Dummy subscript. [VARPRO]

ISYM(600) Integer values of the symbols to be used in plotting data points

(terminal and/or line printer plots). [MAIN, PLTMN, SETUPl, SETUP2,

SPLT, SPLTH, STAT]

ITB Flag to indicate that a table should be printed. [MAIN, STAT]

ITER Number of iterations. [VARPRO]

ITERIN Control variable to determine if a step was retracted. [VARPRO]

ITERM Number of terms in model that depend on a nonlinear parameter.

[MAIN]

ITEST Flag to indicate whether X axis labels are printed in E or F format

(terminal and/or line printer plots). [PLTMN]

ITMAX Maximum number of iterations. [MAIN, VARPRO]

ITP Flag to indicate whether line printer plots should be printed on ter- minal. [MAIN, PR, SETUPl, SETUP2, STAT]

ITRN Flag to indicate that dependent and/or independent variable needs

to be transformed. [MAIN]

145
ITWOM 2*M. [SETUP2]

ITYPE Flag to indicate which type of curve is to be drawn (DISSPLA

plots). [PLOTl] The possible values are:

1 mean line drawn, points sometimes connected;

2 mean line drawn, points never connected;

3 mean line not drawn, points never connected.

IV Number of independent variables per observation. [ADA, CALC, CH,

DERIV, DPA, MAIN, NUMDER, PR, STAT, TRANS, VARPRO]

IVP Number of columns to be printed in table. [PR]

IWT Flag to indicate that weights are calculated in subroutine WT. [MAIN]

IX Subscript used in arranging indexing of IOR. [ORPLOT]

1X1 Plotting location of Ith value of XI vector (terminal and/or line printer

plots). [VMAINH, VPMAIN]

1X2 Plotting location of Ith value of X2 vector (terminal and/or line printer

plots). [VMAINH, VPMAIN]

J Dummy subscript. [ADA, ARRAYD, BACSUB, CH, COV, DPA, MAIN,

NUMDER, 0RFAC1, ORFAC2, ORPLOT, POSTPR, PR, STAT, VARPRO,

VMAINH, VPMAIN]

JJ DUmmy subscript. [ADA, NUMDER]

JM1 J-l (dummy subscript). [COV]

JP1 J+l (dummy subscript). [POSTPR]

JSUB Dummy subscript. [ORFAC1, VARPRO]

K Dummy subscript. [ADA, CALC, DERIV, DPA, MAIN, NUMDER, ORFAC1, ORFAC2, PLTMN, POSTPR, PR, SETUP1, SETUP2, STAT,

VARERR, VARPRO]

KBACK Dummy subscript. [POSTPR]

KG Number of nonlinear parameters. [ADA, CH, DPA, MAIN, NUMDER,

0RFAC1, 0RFAC2, POSTPR, PR, STAT, VARPRO]

KGPK KG+K. [ORFAC2]

KGPKM1 KG+K-1. [0RFAC2]

KGP1 KG+ 1. [0RFAC1, ORFAC2, VARPRO]

KG2 2*KG. [ORFAC21

146
KG23 2*KG+3. [ORFAC1, ORFAC2, VARPRO]

KK Dummy subscript. [ADA, DPA]

KL Dummy subscript. [SETUP2]

KP1 K+l (dummy subscript). [DPA, ORFAC1, 0RFAC2, POSTPR]

KSS(lOl) Integer vector used to determine which plotting symbol is needed

(terminal and/or line printer plots). [PLTMN]

KSUB Dummy subscript. [DPA, VARPRO]

L Dummy subscript. [ARRAYD, COV, DPA]

LABEL Integer value of current row (and Y label) of vertical graphs (terminal and/or line printer plots). [VMAINH, VPMAIN]

LASTC Last column to which Householder reflections should be applied.

[DPA]

LI Dummy subscript. [ARRAYD]

LMSPEC Flag to indicate which set of limits are to be used for multiple,

vertical graphs (terminal and/or line printer plots). [VMAINH, VPLT,

VPLTER, VPLT2, VPMAIN]

LN Flag to indicate whether a label should be printed on the Y axis (terminal and/or line printer plots). [PLTMN]

LVV Integer first dimension in calling program of matrix to be printed.

[ARRAYD]

M Number of observations. [ADA, BAKTRN, CH, DPA, MAIN, ORFAC1, ORPLOT, PLTMN, POSTPR, PR, SETUP1, SETUP2, SPLT,

SPLTH, STAT, TRANS, VARPRO, VMAINH, VPLT, VPLTER, VPLT2,

VPMAIN, WT]

MAX Index corresponding to maximum value of independent variable.

[ORPLOT]

MAXKG Maximum number of nonlinear parameters. [ADA, DERIV,

DPA, MAIN, NUMDER, PR]

MAXN Maximum number of linear parameters. [ADA, MAIN, NUMDER]

MAXNKG MAXN+MAXKG. [MAIN, POSTPR, STAT]

MAXNP2 MAXN+MAXP+2. [ADA, MAIN]

MAXN1 MAXN+1. [ADA, CALC, DERIV, DPA, MAIN, NUMDER, PR,

STAT]

147
MAXP Maximum number of nonzero partials. [CH, MAIN, NUMDER,

PR]

MEANl Plot line location of the XI vector (terminal and/or line printer

plots). [VMAINH, VPMAIN]

MH Highest column number to be printed on a given line. [ARRAYD]

MHMl One less than the single dimension index corresponding to the first element of column MH. [ARRAYD]

MIV Maximum number of independent variables per observation. [MAIN]

ML Lowest column number to be printed on a given line. [ARRAYD]

MLMl One less than the single dimension index corresponding to the first element of column ML. [ARRAYD]

MMAX Maximum number of observations. [ADA, BACSUB, CALC, CH,

COV, DERIV, DPA, MAIN, NUMDER, ORFAC1, ORFAC2, ORPLOT, PLOT2, POSTPR, PR, STAT, TRANS, VARPRO]

MMAX2 2*MMAX. [MAIN]

MMH Single dimension index corresponding to column MH. [ARRAYD]

MML Single dimension index corresponding to column ML. [ARRAYD]

MM1 M-l. [ORPLOT]

MODIT Variable to determine whether information is to be printed at a

given iteration (based on IPRINT). [VARPRO]

MP1 M+l. [SETUP2]

N Number of linear parameters. [ADA, CALC, DPA, MAIN, NUMDER,

ORFAC1, POSTPR, PR, STAT, VARPRO]

NAM(40) Array which reads keywords from saved data file. [MAIN]

NC Flag to skip unnecessary computation. [ADA, MAIN, NUMDER]

NCBLK Common block containing NC. [ADA, MAIN, NUMDER]

NCC Number of columns in matrix to be printed. [ARRAYD]

NCON Number of "constant" functions in A. [DPA]

NCONP1 NCON+1. [ADA, DPA]

NDI(21) Vector of numerical derivative increments. [ADA, CH, MAIN,

NUMDER, PR]

NH Highest row number to be printed on a page. [ARRAYD]

NKG1 N+KG+ 1. [POSTPR]

148
NKG2 N+KG+2. [DPA, MAIN, POSTPR, STAT, VARPRO]

NL Lowest row number to be printed on a page. [ARRAYD, CALC, DERIV]

NM1 N+KG-1. [COV]

NN The number of Y labels to be listed at the left axis (terminal and/or

line printer plots). [PLTMN]

NPH Subscript to identify column number for last function PH. [ADA]

NPK N+K. [ORFAC1]

NPKG N+KG. [COV, MAIN, ORFAC1, POSTPR, STAT]

NPP1 N+P+l. [DPA]

NPP2 N+P+2. [ADA, DPA, VARPRO]

NPTS Number ofpointS to be drawn on graphs (DISSPLA plots). [PLOT1,

PLOT2]

NP1 N+l. [ADA, BACSUB, DERIV, DPA, MAIN, NUMDER, ORFAC1,

POSTPR, PR, STAT, VARPRO]

NP2 N+2. [DPA]

NRR Number of rows in matrix to be printed. [ARRAYD]

NSER Number of vectors to be plotted (terminal and/or line printer plots).

[VMAINH, VPLT, VPLTER, VPLT2, VPMAIN]

NSIZE Flag to indicate graph size (terminal and/or line printer plots).

[PLTMN]

NU Marquardt parameter. [ORFAC2, VARPRO]

NUM Flag to indicate whether all initial values should be printed. [PR]

NUMCOL Number of columns per graph. 101 columns for line printer

plots; 51 columns for terminal plots. [PLTMN]

NUMCP2 NUMCOL+2. [PLTMN]

NUMLAB Number of X labels to be printed at the bottom of graphs

(terminal and/or line printer plots). [PLTMN]

NUMROW Number of rows in graphs (terminal and/or line printer plots).

[PLTMN]

Nl Dimension of triangular matrix used in solving system of equations.

[BACSUB1

149
ONE The plotting symbol "1" (terminal and/or line printer plots). [VMAINH,

VPMAIN]

OUTPUT Unit number to determine where output should be printed.

[POSTPR, VARERR, VARPRO]

P Number of nonvanishing partial derivatives of the PH(J) (i.e., number of ones in INC). [ADA, DPA, MAIN, PR]

PH(ll) Vector containing the values of the functions PH(J) of the model

equation. [ADA, CALC, MAIN, STAT]

PHINP1 Flag to indicate whether PH(NPl) is in the model equation.

[ADA, DPA]

PI 3.141592654 [STAT]

PI4 PI/4. [STAT]

PI8 PI/8. [STAT]

PL(103) Vector of the horizontal print positions for each line of vertical

graphs (terminal and/or line printer plots). [VMAINH, VPMAIN]

PP Variable used in calculating probability of standardized residuals. [STAT]

PRJRES Norm of the projection of the residual onto the range of the

Jacobian of the variable projection functional. [ORFACl, POSTPR,

VARPRO]

R(300) Vector of residuals (Y-ETA). [DPA, MAIN, POSTPR, STAT]

RMS Residual mean square. [STAT]

RNEW Norm of residual after current iteration. [VARPRO]

RNORM Norm of residual. [DPA, POSTPR, VARPRO]

RSD Residual standard error. [STAT]

RSS Residual sum of squares. [STAT]

R2 Coefficient of determination. [STAT]

SAV(300) Vector used to store independent variable values. [ORPLOT]

SAVE Value of A(l,NCON), i.e., first item in last "constant" column of A.

[DPA]

SAVEl Variable used in moving values from one location to another. [POSTPR]

SCALEl Scale of the XI vector for vertical graphs (terminal and/or line printer plots). [VMAINH, VPMAIN]

150
SCALE2 Scale of the X2 vector for vertical graphs (terminal and/or line

printer plots). [VMAINH, VPMAIN]

SIGMA2 Estimated variance, i.e., (RN0RM**2)/(degrees of freedom).

[COV]

SKIP Flag to indicate whether a given iteration should be printed out. [VARPRO]

SM Double precision variable used in computing a sum. [VMAINH, VPMAIN]

SR(300) Vector of standardized residuals. [STAT]

SRV(300) Vector of standardized residuals arranged in ascending order.

[STAT]

STAR The plotting symbol "*" (terminal and/or line printer plots). [VMAINH,

VPMAIN]

SUM Variable used in computing a sum. [COV, STAT]

SWS Sum of the elements of the WEIGHT vector. [STAT]

SYH(300) Vector of probabilities of standardized residuals. [STAT]

SYM(30) Vector of plot symbol assignments (terminal and/or line printer

plots). [PLTMN]

SYM1(9) Vector of integer values used to indicate multiple points at same

location of graphs (terminal and/or line printer plots). [PLTMN]

T(300,6) Array of independent variables. T(I,J) contains the value of the

Ith observation of the JTh independent variable. [ADA, CALC, CH,

DERIV, DPA, MAIN, NUMDER, ORPLOT, PR, SETUP1, SETUP2,

STAT, TRANS, VARPRO]

TEMP(300) Vector of independent variables to be plotted (DISSPLA plots).

[MAIN, STAT]

TITLE(40) Title given to the problem by the user. [MAIN]

TITLED(41) TITLE concatenated with "$" for use as graph lAbel (DISSPLA plots). [MAIN]

TITLE1 First line of graph title (DISSPLA plots). [PLOTl, PLOT2]

TITLE2 Second line of graph title (DISSPLA plots). [PLOT2]

TLABEL(48) Axis label given by user for the independent variable (DISSPLA plots). [MAIN, STAT]

TMAX Maximum value of TEMP vector. [MAIN]

151
TMIN Minimum value of TEMP vector. [MAIN] TPL(600) Vector of values

of the independent variable in the order in which they are to be plotted

(terminal and/or line printer plots). [MAIN, SETUP1, SETUP2]

TR T-ratios. [STAT]

TS Flag to indicate that independent variable is a time series. [MAIN,

PLOT1, PLOT2, STAT]

TWO The plotting symbol "2" (terminal and/or line printer plots). [VMAINH,

VPMAIN]

U(10) Vector containing information on Householder reflections; linear parameters. [DPA, POSTPR]

Ul Factor used in computing Householder reflec- tions. [ORFAC1, 0RFAC2]

VAR Character variable containing "VAR". [MAIN, PR, STAT]

W(300) Vector of nonnegative weights. [DPA, MAIN, POSTPR, VARPRO]

WEIGHT(300) Vector of nonnegative weights, squared. [CH, DPA, MAIN,

PR, STAT, VARPRO, WT]

WRNBEF Flag to control printing certain questions at beginning of "Command Point" only once per execution of the program. [MAIN]

WTVAL User defined value by which the weights for all data points are changed. [CH]

XARAY(300) Vector of observations for the X coordinates (DISSPLA,

terminal and/or line printer plots). [PL0T1, PLTMN, SPLT, SPLTH]

XARAY1(300) Vector of X coordinates for the first curve in multiple

graphs (DISSPLA plots). [PLOT2]

XARAY2(500) Vector of X coordinates for second curve in multiple graphs

(DISSPLA plots). [PLOT2]

XLAB(15) Vector of X axis labels (terminal and/or line printer plots).

[PLTMN]

XLBL(48) X axis label (DISSPLA plots). [PLOT1, PLOT2]

XLIM(9) Vector of the X axis labels for vertical graphs (terminal and/or

line printer plots). [VMAINH, VPMAIN]

XLONG Length in inches of X axis (DISSPLA plots). [PLOTl]

XMAX Upper bound for X axis (DISSPLA plots). [PLOTl, PLOT2]

XMN1 Minimum value of XI. [VMAINH, VPMAIN]

152
XMN2 Minimum value of X2. [VMAINH, VPMAIN]

XMX1 Maximum value of XI. [VMAINH, VPMAIN]

XMX2 Maximum value of X2. [VMAINH, VPMAIN]

XN Minimum X value to be plotted (terminal and/or line printer plots).

[PLTMN]

XORIG Lower bound for X axis (DISSPLA plots). [PLOT1, PL0T2]

XPL(500) Vector of 500 values, equally incremented, of the independent

variable used to plot the predicted Y curve (DISSPLA plots). [MAIN,

STAT]

XWIDTH Interval length of X axis (terminal and/or line printer plots).

[PLTMN]

XX Maximum X value to be plotted (terminal and/or line printer plots).

[PLTMN]

XXLONG X coordinate to determine graph placement on Versatec plotter

(DISSPLA plots). [PLOT1]

XX1(300,18) Matrix product of BM and A(NPKG,NPKG) used in computing standardized residuals. [STAT]

XI (300) Vector of first series of X values to be plotted (terminal and/or

line printer plots). [VMAINH, VPLT, VPLT2, VPMAIN]

X2(300) Vector of second series of X values to be plotted (terminal and/or

line printer plots). [VMAINH, VPLT, VPLT2, VPMAIN]

X1MEAN Mean value of the XI vector. [VMAINH, VPMAIN]

X2MEAN Mean value of the X2 vector. [VMAINH, VPMAIN]

Y(300) Vector of observed values of dependent variable. [CH, DPA, MAIN,

PR, SETUP1, SETUP2, STAT, TRANS, VARPRO, WT]

YARAY(300) Vector of observations for Y coordinate (DISSPLA, terminal and/or line printer plots). [PLOTl, PLTMN, SPLT, SPLTH]

YB Mean of observed values, Y. [STAT]

YC Value of PH(J) using current values of ALF in computing numerical

derivatives. [NUMDER]

YD Value of PH(J) using incremented values of ALF in computing numerical derivatives. [NUMDER]

YFIT(500) Vector of values for predicted Y curve (DISSPLA plots). [PLOT2]

153
YLAB(20) Vector of Y axis labels (terminal and/or line printer plots).

[PLTMN]

YLABEL(48) Axis label given by user for the dependent variable (terminal and/or line printer plots). (MAIN, STAT]

YLBL(48) Y axis label (DISSPLA plots). [PL0T1, PLOT2]

YLONG Length in inches of Y axis (DISSPLA plots). [PLOT1]

YLOWER Lower bound for Y values to be plotted on current line (terminal and/or line printer plots). [PLTMN]

YMAX Upper bound of Y axis (DISSPLA plots). [PLOTl, PLOT2]

YN Minimum Y value to be plotted (terminal and/or line printer plots).

[PLTMN]

YOBS(300) Vector of values for observed Y curve (DISSPLA plots). [PL0T2]

YORIG Lower bound of Y axis (DISSPLA plots). [PLOTl, PLOT2]

YPL(500) Vector of 500 predicted values of the dependent variable, used

to plot the predicted Y curve (DISSPLA plots). [MAIN, SETUPl,

STAT]

YUPPER Upper bound for Y values to be plotted on current line (terminal

and/or line printer plots). [PLTMN]

YWIDTH Interval length of Y axis (terminal and/or line printer plots).

[PLTMN]

YX Maximum Y value to be plotted (terminal and/or line printer plots).

[PLTMN]

YYLONG Y coordinate to determine graph placement on Versatec plotter

(DISSPLA plots). [PLOTl]

Y1MAX Upper bound of first curve for Y axis (DISSPLA plots). [PLOT2]

YIORIG Lower bound of first curve for Y axis (DISSPLA plots). [PLOT2]

Z(10) Solution vector from solving system of equations. [BACSUB]

ZERO(300) Vector of zeros used to draw the zero line on graphs (DISSPLA

plots). [PLOTl]

154
Bibliography

fll Bolstad, John, VARPRO, Computer Science Department, Stanford University, Stanford, California (1974).

Burington, R. S. and D. C. May, Handbook of Probability and Statistics,

McGraw-Hill, New York (1970), pp. 239-240, 243-244.

DISSPLA User's Manual, Integrated Software Systems Corporation,

San Diego, California (1981).

Draper, N. R. and H. Smith, Applied Regression Analysis, Wiley, New

York (1966).

[5l Ehrenberg, A. S. C, Data Reduction, Wiley, London (1975).

[6] Finerty, J. P., The Population Ecology of Cycles in Small Mammals,

Yale Univ. Press, New Haven (1980).

[7] Golub, G. H. and V. Pereyra, The Differentiation of Pseudo-Inverses

and Nonlinear Least Squares Problems whose Variables Separate, SIAM

Jour. Numer. Anal., 10 (1973), pp. 413-432.

[81 Hogben, D., S. T. Peavy and R. N. Varner, OMNITAB II User's Ref- erence Manual, NBS Technical Note 552 (1971), pp. 138-146.

Keith, L. B., Wildlife's Ten-Year Cycle, Univ. of Wisconsin Press,

Madison (1963).

[101 Lawson, C. L. and R. J. Hanson, Solving Least Squares Problems,

Prentice-Hall, Englewood Cliffs, New Jersey (1974).

|llj McKenna, G. B., G. Hadziioannou, P. Lutz, G. Hild, C. Strazielle, C.

Straupe, P. Rempp, and A. J. Kovacs, Dilute Solution Characterization

of Cyclic Polystyrene Molecules and Their Zero-Shear Viscosity in the

Melt, Macromoecules, 20 (1987), pp.498-512.

12] Rail, J. E. and R. E. Funderlic, Interactive VARPRO (INVAR), a Nonlinear Least Squares Program, ORNL/CSD-55, Computer Sciences Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee (1980).

155
[13] Tryon, Peter V. and Janet R. Donaldson, STATLIB, a Library of FORTRAN Subroutines for Statistical Analysis of Experimental Data, Statistical Engineering Laboratory, National Bureau of Standards, Boulder, Colorado (1978).

156

U US GOVERNMENT PRINTING OFFICE: 1987— 181-076/ 61026
FORM NBS-) 14A (REV 1 1-84)

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. NBS/TN-1238

2. Performing Organ. Report No 3. Publication Date July 1987

4. TITLE AND SUBTITLE

An Interactive Nonlinear Least Sauares Program

5. AUTHOR(S)

Claire M. Wolfe, Bert W. Rust, Jeffrey H. Dunn, and Ivor E. Brown

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions)

NATIONAL BUREAU OF STANDARDS

U.S. DEPARTMENT OF COMMERCE

GAITHERSBURG, MD 20899

7. Contract/Grant No.

8. Type of Report & Period Covered

Final

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street, City. State, ZIP)

Same as item 6.

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 significant

bibliography or literature survey, mention it here)

INVAR is an interactive computer code which uses the Stanford variable- separable nonlinear regression program VARPRO to solve nonlinear (and

linear) least squares problems. The variabie-separabie feature of VARPRO

makes it attractive to users with real- world fitting problems because it it- erates only on the parameters which appear nonlinearly in the model. Not

only does this simplify the iteratation. but it also means that the user is not required to supply initial estimates for the parameters which appear

linearly. INNAR implements VARPRO within an environment providing the user with on-line feedback and the opportunity to make changes, transformations, and corrections in real-time. It provides extensive statistical diagnostics and plots of the results. It can provide publication-quality DISSPLA plots either at a graphics terminal during exectution or at a hard copy device

afterwards. The code is written mostly in standard Fortran with only a

few instructions specific to the NBS Cyber S55. With minor changes it can

easily be implemented on other computers.

This report is both a tutorial guide for beginners and a reference manual for experienced users who wish to make changes in the code. It contains

three completly solved example problems. Three appendices contain infor- mation necessary for making changes in the programs.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) interactive computing; least squares; nonlinear regression

13. AVAILABILITY

[jj[] Unl imited

| | For Official Distribution. Do Not Release to NTIS

fx~l Order From Superintendent of Documents, U.S. Government Printing Office, Washington, DC

20402.

^] Order From National Technical Information Service (NTIS), Springfield, VA 22161

14. NO. OF

PRINTED PAGES

160

15. Price

USCOMWDC 85-6004



NBSTechnical Publications

Periodical

Journal of Research—The Journal of Research of the National Bureau of Standards reports NBS research and development 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 time are survey articles on topics closely related to the Bureau's

technical and scientific programs. Issued six times a year.

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 interested 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 Series—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 pro- gram coordinated by NBS under the authority of the National Standard Data Act (Public Law 90-3%).

NOTE: The Journal of Physical and Chemical Reference Data (JPCRD) is published quarterly for NBS by

the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints,

and supplements are available from ACS, 1155 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 perfor- mance 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 re- quirements 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 NBS research and experience, covering areas of interest to the consumer. Easily understandable language and illustrations provide useful background

knowledge for shopping in today's technological marketplace.

Order the above NBS publications from: Superintendent of Documents, Government Printing Office,

Washington, DC 20402. Order the following NBS publications—FIPS 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 Government 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 1 1717 (38 FR 12315, dated May 11, 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 distribution 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

Gaithersburg, MD 20899

Official Business

Penalty for Private Use $300