Numerical Aspects of Some Nonstandard Regression Problems

NOT FOR QUOTATION WITHOUT PERMISSION OF THE AUTHOR

NUMERICAL ASPECTS OF SOME NONSTANDARD REGRESSION PROBLEMS

V. Fedorov A. Vereskov

May 1 9 2 i WP-85-32

Working Papers are interim r e p o r t s on work of t h e International Institute f o r Applied Systems Analysis and have r e c e i v e d only lim- ited review. Views o r opinions e x p r e s s e d herein d o not neces- s a r i l y r e p r e s e n t t h o s e of t h e Institute or of i t s National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 Laxenburg, Austria

ABSTRACT

Regression models a r e extremely popular in d i f f e r e n t areas conjunct with systems analysis. The v a r i e t y of t h e s e models i s immense and now, as a consequence, t h e r e e x i s t many computerized versions of t h e corresponding s t a t i s t i c a l methods. In this p a p e r , we attempt t o unify (from t h e computa- tional viewpoint) a t l e a s t some s t a t i s t i c a l a p p r o a c h e s . We understand t h a t similar attempts have r e p e a t e d l y been made in s t a t i s t i c a l p r a c t i c e ( s e e , f o r instance, BMDP 1983), but none of them c a n b e considered completely suc- cessful. Nevertheless, any new attempt (this one, we hope) gives a more profound and comprehensive understanding of t h e situation and t h e f u t u r e directions of t h e investigations.

CONTENTS

I. INTRODUCTION

11. TRADITIONAL REGRESSION MODEL AND LEAST SQUARE METHOD 111. MAIN LEAST SQUARE (L.S.) ALGORITHMS

IV. MODELS AND ESTIMATORS BASED O N THE LEAST SQUARE METHOD

V a r i a n c e of E r r o r s of O b s e r v a t i o n s Depending upon Unknown P a r a m e t e r s

Estimation of P a r a m e t e r s of Distribution Density Function Multiresponse Case

Regression-Autoregression Models

O b s e r v a t i o n of Deterministic Dynamic System M-Estimates

P r e d i c t o r s S u b j e c t to E r r o r

V. NUMERICAL EXAMPLES

REFERENCES

NUMERICAL ASPECTS OF SOME NONSTANDARD REGRESSION PROBLEMS

V. Fedorov and A . Vereskov

I. INTRODUCTION

The p a p e r gives a s h o r t s u r v e y of models and estimation methods which are closely connected with t h e l e a s t s q u a r e method. The f i r s t two sections a r e devoted t o traditional r e g r e s s i o n models. All o t h e r models and estima- tion methods are considered in t h e t h i r d c h a p t e r . The r e a s o n why t h e main attention w a s paid t o traditional case i s v e r y simple

-

a l l o t h e r models and estimation p r o c e d u r e s u n d e r consideration can b e transformed in some way t o t h i s case. T h e r e f o r e s t a t i s t i c a l ideas and basic p r o p e r t i e s are p r a c t i - cally t h e same f o r all considered s t a t i s t i c a l problems, and more or less detailed analyses p r e s e n t e d in t h e f i r s t two c h a p t e r s allow, f i r s t l y , avoidance of r e p e t i t i o n s a n d , secondly, understanding of some common f e a t u r e s of t h e s e problems.

The final c h a p t e r deals with technical a s p e c t s of 1s-program set (which were p r e p a r e d in VNIISI, Moscow) and some numerical examples a r e p r e s e n t e d in it.

II.

TRADJTIONAL REGRESSION MODEL

AND

LEAST SQUARE METHOD When functional relation between some variable y (usually i t i s called

"response") and variables x

=

( x i ,

- . .

,xk ) (control variables o r predic- t o r s ) must be analyzed, based on empirical d a t a , t h e model

-

y i = v ( x i , a t ) + c i , i = l , n (1)

is traditionally used. In ( I ) , y i is a r e s u l t of a n observation, xi a r e condi- tions of this observation, ~ ( x , 8 ) i s a given function, 8

=

(81,

. . -

^{, 0, ) a r e}

unknown p a r a m e t e r s , s u b s c r i p t " t " stands f o r t r u e value, ci i s a n "error "

of a n observation, "i" i s a number of a n observation. The e r r o r s ci , i

= G ,

c a n d e s c r i b e e i t h e r t h e deviation of function ~ ( z , 0 ) from t h e real function o r "noise" of a n experiment. To justify t h e model (1) i t i s necessary t o make some formal assumptions about t h e i r behavior as well. In traditional c a s e s , i t i s assumed t h a t t h e s e e r r o r s a r e random with z e r o means (E [ci ]

=

^O), independently identically distributed with finite variance (E [cf]

=

b 2 7 2 ( ~ i ), where bZ c a n be unknown and 7 ( x ) is a given function).

Under t h e s e assumptions i t i s reasonable from a s t a t i s t i c a l point of view t o use t h e l e a s t s q u a r e estimation (1.s.e.):

6, =

Arg min TI:(@) ,

8

When ~ ( x ,@) is a l i n e a r function of

8

(e.g., q ( x , B)

= eTf'

( x ) ) t h e minimization problem c a n b e solved explicitly and t h e solution is

n n

where Mn

=

7-2(zi)f'(xi)f'T(xi). Yn

=

7 - 2 ( x i ) ~ i f ' ( ~ i ) . The sum of

i =I i =l

residual s q u a r e s gives t h e estimator of t h e variance:

Z2 =

( n -m)-'v,2(6)

.

This r e s u l t is too well known t o b e mentioned but formula (3) will b e occa- sionally used in what follows. The s t a t i s t i c a l p r o p e r t i e s of 1.s.e. are s u b j e c t t o classical l i n e a r r e g r e s s i o n analysis and c a n b e found p r a c t i c a l l y in any s e r i o u s s t a t i s t i c a l monograph (see, f o r instance, C.R. Rao 1965). F o r t h i s p a p e r i t will b e enough t o mention t h a t

~ ( 6 )

=a2M,I i s a

dispersion ( o r variance-covariance) matrix of

6

^and

- 2 T

d [ ~ ( z , 6 ) ]

=

6 ^f' ( z

)ML'~'

^{( x ) i s a} v a r i a n c e of t h e forecasting:

~ ( z , 6 , )

=

6,Tp(x), at a point z .

In what follows t h e main attention will b e paid t o nonlinear models and fortunately most of t h e p r o p e r t i e s of l i n e a r cases are fulfilled in g e n e r a l situations at l e a s t asymptotically ( s e e Jennrich 1969, Wu 1981).

Let u s s p e l l out t h e principal ones:

-

_{L.s.e. are} (strongly) consistent, i.e., t h e y converge almost s u r e l y t o t r u e value

Bt

when n ⁺

=.

-

The p a r a m e t e r 6' (when 7(xi ) r 1 it i s t h e v a r i a n c e of ci ) con- sistently estimated by

8; =

( n

-

m )-I v,<Qn > .

-

The dispersion matrix of

G ( 6 , -

B t ) i s consistently estimated by

i = l where 'j(z,8)

=

a7j(x, 8)/ a8.

-

The asymptotical distribution (n -, 0 0 ) of

6 ( 6 , - ^{B t )}

i s normal with z e r o mean and v a r i a n c e matrix

n

=

_n lirn _+- n

x

^{f ( x i ,}

^Ot

^{) f '(zit}

^Ot

⁾

i =1

These p r o p e r t i e s a r e fulfilled under r a t h e r mild assumptions which can b e roughly formulated in t h e following way:

-

_Functions 7j(zi, 8) and f ( x i , 8)

=

a7j(ziL, 8) / a@ are smooth enough for any zi in t h e vicinity of

et,

-

The conditions of observations have t o g u a r a n t e e nondegeneracy of t h e estimation problem; for instance, t h e limit function

n

v 2

=

lim n - I

x

^[7j(zi,8)

^-

7 j ( ~ ~ , @ ~ ) ] ~

n +- i = l

h a s t o e x i s t in t h e vicinity of

Bt

and t h e function v2(8) h a s t o have unique minimum f o r 8

= Bt .

Additionally, t h e limit

M(ot) =

lirn n-I n _f (zi,ot)f T(zi,8t)

n -.- _i = I

h a s t o e x i s t and matrix

it(Ot)

h a s t o b e r e g u l a r (rank

M(ot)

= m ) .

The above mentioned p r o p e r t i e s allow us t o use in p r a c t i c e t h e follow- ing approximations:

a n d t h e value

h a s F-distribution with m -q a n d n -m d e g r e e s of f r e e d o m , w h e r e v , z ( ~ , ) r e p r e s e n t s t h e r e s i d u a l sum of s q u a r e s obtained by fitting t h e model with a r e d u c e d number q

<

m of f r e e p a r a m e t e r s .

Minimization problem (2) c a n b e c o n s i d e r e d as a s p e c i f i c c a s e of a non-linear minimization problem a n d a n y a v a i l a b l e g e n e r a l algorithm c a n b e used t o g e t i t s solution

6.

(In what follows index n will usually b e omitted).

However, as w a s r e p e a t e d l y discussed in s t a t i s t i c a l l i t e r a t u r e ( s e e , f o r i n s t a n c e , Chambers 1 9 7 3 , Fedorov e n d Uspensky 1 9 7 5 , J e n n r i c h a n d Ralston 1978), d u e t o t h e s p e c i f i c s t r u c t u r e of v 2 ( 8 ) s o m e algorithms are much more e f f e c t i v e t h a n g e n e r a l ones. Moreover, o n e should k e e p in mind t h a t b e s i d e s e s t i m a t e s

6

f o r a n y s t a t i s t i c a l analysis, i t is n e c e s s a r y t o c a l c u l a t e m a t r i c e s ~ ( 6 ) a n d function d[7j(z, 8 n ) ] f o r a p r e s c r i b e d set zl,

.

_sXnp.

Most of t h e a l g o r i t h m s used in s t a t i s t i c a l p r a c t i c e c a n b e r e p r e s e n t e d i n t h e form

Qs

=

Qs + P, Hs Ys ( 6 )

w h e r e s is t h e number of a n i t e r a t i o n , p, i s t h e length of a s t e p , Hs i s some m a t r i x (positively semidef i n i t e in most c a s e s ) ,

If H, E I t h e n (6) i s a n algorithm of t h e g r a d i e n t t y p e . When

H,

=

(M, + Y,A,)-' , w h e r e y,

>

^{0 ,} A, r 0 , a n d

o n e d e a l s with algorithms of t h e Gauss-Newton t y p e . The r e g u l a r i z a t o r s y, a n d A, are usually non-zero f o r ill-conditioned minimization problem (2).

Algorithms of Gauss-Newton t y p e are p r o b a b l y most p o p u l a r in s t a t i s t i - c a l s o f t w a r e . They c o n v e r g e much f a s t e r t h a n algorithms of t h e g r a d i e n t t y p e a n d d e l i v e r some additional information ( f o r i n s t a n c e m a t r i x M;') which i s n e c e s s a r y f o r s t a t i s t i c a l analysis. But t h e y demand calculation of d e r i v a t i v e s f ( z i , 8,) at e v e r y s t e p (by t h e way, a n y algorithm (6) demands it). Unfortunately, in modern e m p i r i c a l r e s e a r c h , t h e calculation of func- tions ~ ( z ,@) and f ( z , 0) h a p p e n s t o b e a s e r i o u s numerical problem ( f o r i n s t a n c e , ~ ( z , 8) c a n b e d e s c r i b e d by some system of d i f f e r e n t i a l equations) or v e r y tedious f o r a p r o g r a m m e r . Of c o u r s e , a s t r a i g h t f o r w a r d u s e of fin- i t e d i f f e r e n c e approximation

where e, i s a v e c t o r with z e r o components e x c e p t t h e a-th o n e , c a n s a v e from tedious programming b u t i t i s n o t r e a s o n a b l e f r o m numerical point of view.

In g e n e r a l nonlinear programming, t h e r e is a v a r i e t y of d i f f e r e n t methods f o r avoiding calculation of d e r i v a t i v e s ( s e e , f o r i n s t a n c e , Chambers

1973, Fedorov and Uspensky 19'75, Himmelblau 19'73). The most efficient among them are based on d i f f e r e n t q u a d r a t i c approximations of v 2 ( @ ) . But a l l of them were d e f e a t e d in s t a t i s t i c a l applications by methods which are based o n l i n e a r approximation of ~ ( z , 8) using t h e "history" of i t e r a t i v e p r o c e d u r e in conjunction with t h e idea of Gauss-Newton method.

In t h e s e methods (presumably, f i r s t suggested by Peckham ( 1 9 7 0 ) ) , at e v e r y s t e p of i t e r a t i v e p r o c e d u r e t h e following minimization problem must b e solved.

s'=s -1

f i s

=

Arg min

C

o S s r [ v ( z i , 8 , . ) - 7 7 ( ~ ~ , @ , - ~ )

- ~ T ( @ s * -

@s-1)I2 ¹ fi ^,'=(I

w h e r e u s , , are some weights. In t h e simplest case ( s e e Ralston and

- - . . .

Jennrich 19'78) u s , ,

-

^{u s , , -2}

- =

1 a n d a l l rest weights equal to z e r o . Similarly t o (2) a n d (3). o n e c a n g e t t h a t

w h e r e

I t i s r e a s o n a b l e t o s u g g e s t t h a t approximately

P

( x i *

os

-1)

=

.Pis ^,

N o w o n e c a n u s e i t e r a t i v e p r o c e d u r e (6) and ( 7 ) with

A f t e r some e l e m e n t a r y t r a n s f o r m a t i o n , t h i s p r o c e d u r e c a n b e p r e s e n t e d in a more convenient f o r m :

n

where

& = z

y-'(xi )Uis

~ f T s

i =1

The p r o p e r t i e s of i t e r a t i v e p r o c e d u r e (9) w e r e studied by Vereskov ( 1 9 8 1 ) and i t should b e pointed o u t t h a t they depend upon not only s t r u c t u r e

of ~ ( x , 0 ) and sequence x l,

.

, x , but also upon sample c l ,

- -

^{, E , . In}

o t h e r words, a l l a s s e r t i o n s about convergency, for instance have a proba- bilistic c h a r a c t e r and t h e r e e x i s t with non-zero probability s o m e samples when t h e p r o c e d u r e does not converge in s p i t e of "good quality" of ~ ( x ,@)

The i t e r a t i v e p r o c e d u r e ( 9 ) i s t h e basic p a r t of a l l algorithms dis- cussed below and realized in program LsO.

IV.

MODELS

AND

ESTIMATORS BASED ON THE LEAST SQUBRE METHOD

V a r i a n c e of E r r o r s o f O b s e r v a t i o n s D e p e n d i n g upon Unknown P a r a m e t e r s

In t h e traditional c a s e , i t is assumed (see Ch. 2 ) t h a t t h e function y 2 ( x i ) is known. In t h e more g e n e r a l case, i t i s n a t u r a l t o assume t h a t vari- a n c e of random values ci depends upon unknown p a r a m e t e r s :

If a l l p a r a m e t e r s in ( 1 0 ) coincide with some p a r a m e t e r s of r e s p o n s e func- tion ~ ( x ,@) t h e n t h e estimator defined by t h e i t e r a t i v e p r o c e d u r e :

is a c c e p t a b l e f r o m a s t a t i s t i c a l point of view ( s e e Fedorov 19'74). To solve ( 1 2 ) one can use i t e r a t i v e p r o c e d u r e (9). The corresponding algorithm is

r e a l i z e d in p r o g r a m Lsl.

If function 7 2 ( z , 8 ) contains p a r a m e t e r s which are n o t involved in

~ ( z ,

a ) ,

t h e n a m o r e complicated p r o c e d u r e i s needed to g e t t h e e s t i m a t o r of unknown p a r a m e t e r s ( s e e Malyutov 1982):

n

Bq

=

Arg min

x

^{{7-'(zi, Oq} - l ) [ ~ i

-

l ) ( z i . @)I' +

@ i = 1

The p r o p e r t i e s of ( 1 3 ) f o r normally d i s t r i b u t e d ci are studied by Malyutov (1982). In o t h e r cases, e s t i m a t o r ( 1 3 ) i s s t i l l c o n s i s t e n t a n d asymptotically normally d i s t r i b u t e d u n d e r some r e a s o n a b l y mild conditions.

Naturally, o n e c a n u s e p r o g r a m LsO to solve ( 1 4 ) . The c o m p u t e r r e a l i - zation of ( 1 3 ) a n d ( 1 4 ) i s included in t h e system u n d e r t h e name Ls2.

Estimation o f Parameters of Distribution Density Function

L e t z b e a random value with density p ( z , @ ) , z E X , 8 ^E

fl

C R m . The i - t h e x p e r i m e n t c o n s i s t s of r i o b s e r v a t i o n s of numbers of cases nij when z E

q j

c X ,

Xij ^u q1 =

0, j f L

.

The method of estimation of @ closely con- n e c t e d with t h e 1.s. method was suggested by Vereskov a n d Pshennikov (1983). They a l s o discussed i t s links with t r a d i t i o n a l a p p r o a c h e s . The c o r r e s p o n d i n g numerical p r o c e d u r e h a s t h e following form:

n _{I t} Ni [ u i j

-

p U ( @ ) l Z Bq

=

Arg min

x x

@ i = 1 j = 1 ~ i j ( @ q - l )

T i T i

w h e r e Ni

=

nil , p i j (0)

=

/ P ( 2 . @ ) d l /

x

^{/ P ( 2 .} @ ) d z l uij

=

^nij / N i .

i =l

xv

j = q j

In t h e simplest case o n e c a n use approximation p i j ( 8 ) ' p ( x i j , @ ) b i j , d i j

=

I d x . In g e n e r a l c a s e , i n t e g r a l s a r e c a l c u l a t e d numerically. The

x,

algorithm i s r e a l i z e d in p r o g r a m ls3.

Multiresponse Case

If in model ( 1 ) t h e r e s p o n s e i s a v e c t o r , y i E

R'

, a n d E [ c i ciT]

=

Di , w h e r e Di are given f o r a11 i

= l,n

t h e n t h e solution of t h e minimization problem

6 =

Arg min

x

n ^{[ p i}

^-

r l ( x i , O ) I ~ D ~ - ' [ Y ~

-

r l ( x i ,

@ ) I

^{( 1 5 )}

i = 1

c a n b e used as t h e e s t i m a t e of 8. I t is c l e a r t h a t ( 1 5 ) h a s t h e same s t r u c - t u r e as (2). The c o r r e s p o n d i n g numerical p r o c e d u r e i s contained in p r o - gram L s 4.

If t h e d i s p e r s i o n m a t r i c e s Di are unknown b u t t h e r e is t h e p r i o r infor- mation t h a t

E [ E ~ E [ ] = D ( O t )

.

r a n k D ( e t )

=

¹ ,

t h e n t h e e s t i m a t o r of 8 c a n b e d e f i n e d by t h e i t e r a t i v e p r o c e d u r e (Fedorov 1977, Phillips 1976)

Under some mild condition

6

is s t r o n g l y c o n s i s t e n t and asymptotically r a n - dom values

6 ( 6 -

O t ) are normally d i s t r i b u t e d . Moreover,

is a consistent estimator of D ( B t ) .

P r o c e d u r e ( 1 6 ) i s realized in program Ls5.

Regression-Autoregression Models

In t h i s c a s e

pi

=

^{7 ) ( q} , p i -1, ' ' ' ,Yi - k , 0 ) + ei

.

( 1 7 )

Formally, v a r i a b l e s p i

. .

_{, y i -k} can b e joined t o t h e set of independent variables xi ( s e e , f o r instance, Anderson 1971) and t h e values

can be used as estimates for 0 .

The problem ( 1 8 ) practically coincides with ( 2 ) but f o r convenience only set ( y i ,xi

)In

^{i s used} as an input in Ls ₆

Observation of Deterministic Dynamic System

In systems analysis v e r y often a r e s p o n s e function ~ ( x , 0 ) can b e described by a system of o r d i n a r y differential equation

This specific c a s e of r e g r e s s i o n model ( 2 ) can b e t r e a t e d with t h e help of Ls 7 .

When a r e s e a r c h e r s u s p e c t s t h a t between random v a l u e s ci c a n b e out- lined t h e so-called r o b u s t estimation methods are recommended. One of t h e most p o p u l a r methods i s M-method and c o r r e s p o n d i n g e s t i m a t o r s are defined by t h e minimization problem

w h e r e p(

I

z

/

⁾ i s usually some monotonous n o n d e c r e a s i n g function of ( z j The solution of ( 2 0 ) u n d e r some mild conditions ( s e e , f o r i n s t a n c e , Mudrov a n d Kushko 1976) c a n b e found with t h e h e l p of t h e i t e r a t i v e p r o - c e d u r e :

O .. =

lim Bq ,

q ^+-

P [

I

y i

-

T ( x ~ , Qp -1

I I

Oq

=

Arg min [ p i

-

^77(xi

, @ ) 1 2

i'l

-

^{7 ) ( q}

,a,,

w h e r e t h e a u x i l i a r y minimization problem i s t h e 1.s. problem. Usually t h e stabilization of Bq h a p p e n s a f t e r 3-4 i t e r a t i o n s . P r o c e d u r e ( 2 1 ) i s r e a l i z e d in Ls8 comprehensive discussion of s t a t i s t i c a l p r o p e r t i e s of ( 2 0 ) c a n b e found, f o r i n s t a n c e , in E r s h o v (1978), H u b e r ( 1 9 7 2 ) .

Predictors Subject t o Error

If a r e s e a r c h e r wished t o o b s e r v e h i s system u n d e r conditions xi but b e c a u s e of some random impact t h e y happened t o b e ui

=

xi

+

h i , t h e n , instead of ( 1 ) o n e d e a l s with t h e model

y i

=

v ( x i

+

hi ,@)

+

r i , i

= -

1 , n , ( 2 2 ) w h e r e a l l ci and hi are independent a n d E [ $ ]

=

7 2 ( z i ) , E [ h i h:]

=

D ( x i ).

For model (22) t h e following estimation can b e used (see Fedorov 1974):

6 =

lim Oq ,

q +-

TL

Oq

=

Arg min X(zib Oq - l ) [ ~ f

-

^%zit

@)I'

0 f = l

where

The p r o c e d u r e (23) i s realized in l s 9 .

In conclusion w e emphasize once more t h a t a l l algorithms described in t h i s section are based on IsO.

V.

NUMERICAL

EXAMPLES

To illuminate t h e possibilities of t h e considered s e t of programs t h r e e simple regression problems will b e considered. The d a t a a r e borrowed from t h e p a p e r by C. Marchetti, (1983) (see printout I ) , and a r e extended f o r a f e w additional y e a r s . They d e s c r i b e t h e c a r population in Italy. In t h e cited p a p e r i t w a s suggested t o use a r e g r e s s i o n model with logistic response function and with additive uncorrelated random errors ( t h e l a t t e r statement is not e x p r e s s e d explicitly but it follows from t h e context):

where 8

=

( o l , Q 2 , a r e unknown p a r a m e t e r s , zf stands f o r time.

Formula (24) does not d e s c r i b e completely t h e r e g r e s s i o n model. I t i s still necessary t o clarify t h e assumed p r o p e r t i e s of e r r o r s e f .

T h r e e possible v a r i a n t s will b e considered h e r e : 1) ^{Variance E[E:]} is independent of zi and constant.

2) Variance E [E:] is equal t o 62.r12(zf, @), where 6' h a s t o b e a l s o estimated ( d E [ c f 2 ] / .rl(zf. @)

=

const. "relative e r r o r " is con- s t a n t ) .

3) ~ [ r : ]

=

62.r1(zi,@).

In t h e f i r s t case, program Is0 was used. The r e s u l t s are evident from prin- tout 1.

I t is interesting t o stress two facts. Firstly, t h e residuals (see column 'Y-F" and "NO.RESM and comments in t h e printout) have a tendency t o i n c r e a s e . Secondly, t h e i r signs a r e definitely not randomly distributed.

T h e r e f o r e , o n e may s u s p e c t t h a t t h e more complicated second c a s e is c l o s e r t o r e a l i t y and moreover t h e e r r o r s are c o r r e l a t e d or r e s p o n s e func- tion does not r e f l e c t r e a l i t y . W e are not concerned h e r e with a c c u r a t e sta- tistical analysis of t h e problem b u t only with illumination of how t h e software is working and t h e r e f o r e w e r e s t r i c t o u r s e l v e s only by s t r u g g l e with nonhomogeneity of e r r o r s using t h e hypothesis of cases 2 and 3.

The r e s u l t s are on printout 2. A s initial values f o r estimated parame- t e r s , t h e estimates f r o m t h e previous case were t a k e n . The r e l a t i v e discrepancy between estimates happened t o b e more t h a n 5%. Of c o u r s e , this is not t o o much b u t is s e v e r a l times more than t h e i r s t a n d a r d e r r o r s . T h e r e f o r e one c a n assert t h a t t h e correponding two models give signifi- cantly (in a s t a t i s t i c a l s e n s e ) d i f f e r e n t r e s u l t s .

Unfortunately f o r t h e second model, t h e r e s i d u a l s have a n i n v e r s e ten- dency: t h e y d e c r e a s e in a v e r a g e . Recollecting t h a t in growth p r o c e s s e s t h e

o b s e r v e d value v e r y o f t e n d i s t r i b u t e d a c c o r d i n g t o Poisson's l a w , t h e t h i r d v e r s i o n with

=

v ( z i , 8) was analyzed. The f i n a l r e s u l t s are in p r i n t o u t 3. I t i s c l e a r t h a t now r e s i d u a l s h a v e no t e n d e n c y t o i n c r e a s e o r d e c r e a s e systematically b u t t h e i r signs a p p e a r in long s e r i e s .

I t seems t h a t a l l t e c h n i c a l d e t a i l s on programming are c l e a r f r o m t h e p r i n t o u t s a n d marginal comments. Usually t h e c o n t e n t of input information a r e defined by questions which a p p e a r on t h e s c r e e n a f t e r a p p l i c a t i o n t o a p r o g r a m f r o m t h e c o n s i d e r e d set. S u b p r o g r a m s f o r r e s p o n s e functions a n d weight functions should b e l o c a t e d in a u x i l i a r y f i l e s "resp" a n d "weight,"

c o r r e s p o n d i n g l y . F o r o t h e r p r o g r a m s i t becomes n e c e s s a r y t o u s e some additional a u x i l i a r y files. F o r i n s t a n c e , p r o g r a m ls3 u s e s f i l e

'?>EN"

f o r a density function p ( x , B ) , p r o g r a m l s 6 u s e s f i l e "AVT" f o r a u t o r e g r e s s i o n function, p r o g r a m l s 7 u s e s f i l e "DIFUR" f o r \k(z,B) ( s e e (19)) a n d s o on.

More detailed information c a n b e obtained f r o m IIASA's Computer S e r v i c e s .

s e e I

Note N o . 1

I D A T A : Y , X , W

N a l c (1) ( 1 )

i

1 0.3428e+00 0 . e+0a 0 .1000e+Ol

I 2 a .6130e+00 0 -3000eK11 0.10OOeM1

I 3 0 .6910e+00 0 . 4 0 0 0 e M l 0.1000eK11

I 4 0.8610e+00 0.5000e+01 0 . 1 0 0 0 e M l

I 5 0.1031e+01 0.6a110eMl Q . 1 0 0 0 e t 0 1

I 6 0.1231e+01 0.7000e+01 O.l0OQe+01

I 7 0 - 1 3 9 3 e M 1 0 .B300et01 O .li130e+01

I 8 0 . 1 6 5 9 e 4 1 0 . 9 0 0 Q e + a l a .lOr33e+31

I 9 0.197GeMl 0 . 1 0 0 0 e + 0 2 0 .lOa(le+Ol

I 10 0 .2449e+01 0 .1100e+(32 0 .ldOOe+0l

I 1 1 0 .3030e+01 0 . 1 2 0 0 e M 2 (3 . 1 0 3 0 e + a l

1 12 0 . 3 9 1 3 e M l 0.1300e+02 0 . 1 0 0 0 e M l

1 1 3 0.4675e+01 0.14O(le+02 0 . 1 8 0 0 e 4 l

1 14 0 . 5 4 7 3 e M l 0.1513Oe+02 O .lr330e+Ql

1 1 5 Q.G357e+01 0.1600e-2 O . l B J 0 e 4 1

1 16 0.7295e+01 0.1700eK12 0.1090e+01

1 17 0.8266e-1 8 . 1 0 0 0 e + 0 2 0 . 1 0 0 3 e + 0 1

1 18 0 . 9 1 7 4 e M l 0.1900e+02 0.1@30e+Ol

1 1 9 0.1Q19e+02 0 . 2 0 0 0 e M 2 0 . 1 0 0 0 e + 0 l

1 ²⁰ 0 .1129e+02 0 . 2 1 0 0 e M 2 0 .lllllUe+Ol

1 21 0 . 1 2 4 8 e + 8 2 0 . 2 2 0 0 e M 2 O.l000e+01

1 22 0 . 1 3 4 3 e M 2 0.2300e+02 Q .1000e+01

1 ²³ 0 . 1 4 3 0 e M 2 0.2400e+02 0 . 1 0 0 0 e M 1

1 24 0 .1506e+02 0 .250ae+02 0 .10a0eK11

1 25 0.1593e+02 0.2600e+02 0 . 1 0 a 0 e + 0 l

1 26 0.1647e+02 0.2708e+02 (I - 1 Q 0 0 e M l

1 27 0 . 1 6 2 4 e M 2 0.28(lOe+02 O . 1 0 3 Q e + a l

1 28 0 . 1 7 1 3 e M 2 0 . 2 9 0 0 e M 2 0 . 1 0 0 0 e M l

1 29 0 . 1 7 0 2 e M 2 0.3000e+02 0.1000eii31

1 30 0.1778eK12 0 . 3 1 0 0 e M 2 0 .1000eK11

1 31 0 .1845e+02 0 .3200e+02 O .1Q0Oe+01

I

( 2 ) I NUMBER O F P A M T E R S 3

I NUNBER O F V A R I A B L E S 1

I NUEIBER O F C A S E S 3 1

I I N I T I A L PARAMETERS -0.1000e+32 Q .500Oe+00 0 . 2 5 0 0 e M 2

I

D a l a should b e s a v e d In Lhe a u x i l i a r y f i l e "enl.dn1a". T h e f l r s l column i s o h s e r v a l i o n s ( d e p e n d e n l v a r l a b l e o r r e s p o n s e ) . T h e n e x l c o l u m n s c o n t a i n t h e v a l u e s of p r e d i c t o r s ( i n d e p e n d e n t v a r i a b l e s , c o n l r o l s ) . T h e l a s l column d e s c r i b e s "weiphls" of o b s o r v a l l o n s : W[ =

oc2

in Lhe s i m p l e s 1 c a s e . When only r a l i o s b e l w e e n v a r i a n c e s a r e known, l h e n W[ = & L 7 - 2 ( 2 [ ) , w h e r e &L

will b e e s t i m a t e d a n y Initially a n d r e a s o n a b l e p o s l l i v e n u m b e r c a n b e u s e d In Inpul s l a l e m e n l ( d e l a i l s s e e In Lhc main Lev1 or In s u b s e q ~ i e n l columns). If o n e wish t o g e l f o r e c a s t i n e a t s o m e prescribed p o l n l s whlch d o no1 b e l o n g Lo

"daLa" t h e s e p o i n l s may b e l n l r o d u c e d wllh "small" w e l g h l s (1.e.. lo4).

N a l e (2)

Number o f e s l i m a l e d p a m m e l e r s , n u m b e r of l n d e p e n d e n l variables. n u m b e r of o h s e r v a l i o n s , l n i l l a l v a l u e of e s l l m a l e d p a r a m e l e r s .

TARA?IE'I'ER ERROR = 0

I.IAiII.IUM NU!,:RER O F I T E M ' T I O I J S 5 0 t J U I : X R O F F R E E PhF.hFlETERS 3 l f 4 E I R t i U b I E E R S 1 2 3

DELTA 1 N I T I T A 11= O . l O e + O Q , DELTA LAST H 1 = 0 . 1 0 e 4 1 L I K I T FOR P I V O l ' I t I G : T O k O . 1 0 e - 0 9

NUMBER O F D I V I S I O N S : K 1 = 2

NUI.IEER O F D I V I S I O 1 : S F O R RANDOM VECTOR: K 2 = 2 CONSTANTS F O R CONVERGENCE C R I T E R I O N : L 1 = 2 , L 2 = 2

REDUCTION O F RES.SUFi: REDS= 1 . 3 3 0 d ! I e - 1 U 8

I T E R . ADD. R E S - S U M PARAMETERS

kote (3.1)

2

If only t h e r a t i o s b c t w e e n v a r i a n c e s of o b s e r v a l i o n s a r c known.

t h e n 0 s h o u l d b e u s o d . If all v a r i a n c e s ( o r weip,hLs) a r c known.

t h e n 1 s h o u l d b e p u t In.

5

w

I

Note (3.2)

Somelimcs i t i s u s e f u l to fix s o m e of t h e e s t i m a t e d p a r a m e t e r s .

8

T h e n t h e i r - n u m b e r s s h o u l d n o t a p p e a r h e r e .

C

^c:

Note (3.3)

2

D e f i n e t h e s i z e of o p c r a b l l f t y r e g i o n of l i n e a r a p p r o x i m a t i o n ( s e e p a g e 7 ) .

Note (3.4)

I J p p e r bound f o r a n u m b e r of s l e p - l e n g t h divisions.

Note (3.5)

U p p e r bound f o r a n u m b e r of h a l f i n g s of t h e a d d i t i o n a l r a n d o m v e c t o r l o n g l h . Usually moving in s o m e randomly c h o s e n d l r e c t i o n h e l p s t o a v o i d s l n g u l a r i l y of a p p r o x i m a l l o n (9).

Note (3.6)

t h e n u m b e r of I t e r a t i o n s when t h e I n e q u a l i t y v:

2 < ld.' must b e fulfilled (v,? =

5

^{Wiryt - q ( z i .@,)]2)}

t =1 a n d t h e p r o g r a m wlll b e t e r m l n a t e d .

Note (3.7)

Additional c h o l c e of t h e t e r m i n a t i n g of t h e p r o g r a m . T h e p r w p r a m wlll b e t e r m f n a t e d If v i / v : >REDS.

Printout 1

-

continued

RESIDUAL SUM / ( N - M )

a . 3 8 6 7 0 2 e - 0 1

STANDhRD DEVIATIONS OF PAWLMETER EST1t:ATES 0 . 6 2 6 d - 0 1 0 . 4 1 4 d - 0 2 0 - 1 6 3 d + 0 0 ESTIPUITE OF COVRRlANCE MATRIX

0 . 3 9 1 d - 0 2

- . 2 4 7 d - 0 3 8 . 1 7 2 d - 0 4

0 . 6 2 5 d - 0 2 - . 5 3 5 d - 0 3 0 . 2 6 5 d - 0 1 ESTIFVITE OF CORRELATlON MATRIX 0 . 1 3 8 d M 1

- . 9 5 3 d + Q 0 0 . 1 0 0 d + 0 1

8 . 6 1 4 d c 3 0 - . 7 9 3 d + a 0 0 . 1 8 3 d + B 1

H l S T O G W OF KESlDUhLS Y-F

NoLe ( 5 ) ~.

E s l i m a t r s of unknown p a r a m e t e r s and Lht:ir c o v a r i a n c e and c o r r e l a l i o n m a t r i c e s . I1 i s uselul Lo remember LhaL residual surn/(N-M) = 2 , s r e a l s o

NoLC? ( 1 ) .

8

_a

NoLe (6.1)

Addit.ional oritputs which c a n htr useful In t h e modcl Leslinp and f o r c o a s l i n p

T E S T 61: 1F:CCFCt:DEt:CE O F K E S l D L ' P L S

IJUWBER O F P O I N T D I S T N J S E Y- F

T A B L E I

Note (6.2)

5

-- 7 - T N w

D I S T A N C E = d ( z , - z) ( r , - I ) , G = N-' ^{z , .} This InformaLion Is usrful I

I = I

in Lhe c a s e or mulLidimensiona1 z .

8

u

E

&

P R O G R A M L S 1

DATA: Y,X

NUMBER O F PARPMETERS 3

NUMBER O F V A R I A B L E S 1 NUMBER O F C A S E S 31

I N I T I A L PARAElETERS - 0 . 4 3 0 0 e + 0 1 0 . 2 2 5 0 e + a 0 0 . 2 3 0 0 e + ( 1 2

I

l N T E R N A L CONSTANTS N

I

Kumber n l Iterations wlthln basic t s O . Usually It s h o ~ l d h c equal to 3-5 ^P _I

PARRFlETER ERROR = 0

PUS IEWM NUMBERS O F I T E R A T I O N S

NUMBER O F F R E E PARAMETERS 3

T I l E I R NUMBERS 1 2 3

5a

a

^{Number of} iterations d s s c r i h e d by (11). (12). .ape 8

DELTA I N I T I A L H = 0 . 3 0 e - D l , DELTA L A S T H 1 = 0 . 5 0 e - 4 3 2 L I M I T FOR P I V O T I N G : T O G - 0 . l O e - 0 9

NUMBER O F D I V I S I O N S : K 1 = 2

NUMBER O F D I V I S I O N S FOR RANDOM VECTOR: K 2 = 2 CONSTANTS FOR CONVERGENCE C R I T E R I O N : L 1 = 2 , L 2 = 3 I T E R . ADD. R E S . S U M P A M E T E R S

Tk1E R E S I D U A L SUM W*(Y-F)**2 THE RESIDUAL SUM /(N-M)

0.784476e-01 0.280170e-02

ESTIMATE O F PARAMETERS

-0 . 4 1 6 3 3 6 e M l 0 .207376e+00 0 .202167e+02

STANDARD D E V I A T I O N S OF PARAMETER ESTIMATES 0 - 3 1 2 d 4 l 0 - 3 1 2 8 - 0 2 0 .550d+00

ESTIMATE O F COVARIANCE MATRIX 0.9768-0 3

-.23Od-04 0.975d-55

- .826d-02 - .l15d-02 0 . 3 0 2 d M 0

S T A . D E V . F Nl

ESTIMATE O F CORRELATION MATRIX 0.100d+01

- - 2 3 5 d i 0 0 0 . 1 0 0 d M 1

- .481d+00 - ^.668d+30 0 . 1 0 0 d M l

THE RESIDUAL SUN W f ( Y - F ) * * 2 THE RESIDUAL SUM / ( N - M ) 0 . 1 7 9 2 5 8 e + 0 0 0 - 6 4 0 2 0 6 e - 0 2 E S T 1 MATE OF PARAMETERS

- 0 . 4 2 4 1 2 5 e + 0 1 0 . 2 1 9 3 9 4 e 4 0 0 . 1 9 2 3 2 1 e 4 2 STANDARD DEVIATIONS OF PARAMETER ESTIMATES

0 - 3 7 8 d - 0 1 0 . 3 0 9 d - 0 2 0 . 2 2 2 d 4 8 ESTIMATE O F COVARIANCE MATRIX

0 . 1 4 3 d - 0 2

- . 9 8 9 d - 0 4 0 . 9 5 6 d - 5 5

0 - 2 6 2 d - 3 2 - . 4 9 0 d 4 3 0 - 4 9 4 8 - 0 1

1. RES

ESTIMATE OF CORRELATION FlATRM 0 . 1 0 0 d W l

- - 8 4 7 d W 0 0 . 1 0 0 d + 0 1

0 . 3 1 2 d U 3 a - . 7 1 4 d + 0 0 0 . 1 0 0 d + 0 1

REFERENCES

Anderson, T.W. 19'71. The S t a t i s t i c a l A n a l y s i s of Time S e r i e s . N e w York:

Wiley.

BMDP S t a t i s t i c a l Software. 1983. W.J. Dixon, e d i t o r . University of Califor- nia P r e s s .

Chambers, J.M. 19'73. Fitting nonlinear models: numerical techniques.

B i o m e t r i k a , 60:3-18.

Ershov, A.A. 19'78. Robust method of estimation. A u t o m a t i k a a n d Telemechanics, 8:25-50.

Fedorov, V.V. 19'74. Regression problems with controllable v a r i a b l e s sub- ject t o e r r o r . B i o m e t r i k a , 61:49-56.

Fedorov, V.V. and A.B. Uspensky. 19'75. Numerical a s p e c t s of l e a s t s q u a r e s method in design and analysis of r e g r e s s i o n experiments. Moscow State University (in Russian).

Fedorov, V.V. 19'7'7. Estimation of r e g r e s s i o n p a r a m e t e r s in multiresponse c a s e . In R e g r e s s i o n E x p e r i m e n t s . Moscow: Moscow University (in Russian).

J e n n r i c h , R.I. 1969. Asymptotic p r o p e r t i e s of nonlinear least s q u a r e s esti- mator. A n n . Math. S t a t i s t . 1, 40:633-643.

J e n n r i c h , R.I. and M.L. Ralston. 19'78. Fitting nonlinear models t o d a t a . Technical R e p o r t No. 46. Los Angeles: University of California.

Huber, P.J. 1972. Robust statistics: a review. A n n . Math. S t a t i s t . , 43:1041-1067.

Malyutov,

M

.B

.

1982. Asymptotical p r o p e r t i e s and applications of IRDEIN estimators of p a r a m e t e r s in g e n e r a l r e g r e s s i o n model. In S t o c h a s t i c Process a n d AppLications. pp. 144-158. Moscow (in Russian).

Marchetti, C. 1983. The automobile in a system context: t h e p a s t 8 0 y e a r s and t h e next 20 y e a r s . RR-83-18. Laxenburg, Austria: International Institute for Applied Systems Analysis.

Mudrov, V.I. and V.L. Kushko. 1976. Methods of A n a l y s i s of O b s e r v a t i o n s . Moscow: Soviet Radio (in Russian).

Peckham, G . 1970. A new method f o r minimizing a sum of s q u a r e s without calculating gradients. C o m p u t e r J., 13:418-420.

Phillips, P.C.B. 1976. The i t e r a t e d minimum distance estimator and t h e quasi-maximum likelihood estimator. Econometrics, 44:449-460.

Ralston, M.L. and R.I. Jennrich. 1978. DUD, a derivative-free algorithm f o r nonlinear least s q u a r e s . Technometrics, 20:7-14.

Rao, C.R. 1965. L i n e a r S t a t i s t i c a l Methods a n d T h e i r A p p l i c a t i o n s . N e w York: Wiley.

Vereskov, A.I. 1981. Convergence of some d e r i v a t i v e f r e e algorithms. In Methods of Complex S y s t e m s S t u d y i n g . pp. 58-62. Moscow: VNIISI (in Russian).

Vereskov, A.I. and A.S. Pshennikov. 1983. Estimation of a g g r e g a t e d param- eters of work duration. In D y n a m i c s of Non-Homogeneous S y s t e m s . Moscow: VNIISI (in Russian).

Wu, C.F. 1981. Asymptotic t h e o r y of nonlinear l e a s t s q u a r e s estimations.

Ann. S t a t . , 9:501-513.