NUMAL Section 5.2.1.3.1

BEGIN SECTION : 5.2.1.3.1 (February, 1979)

AUTHOR : B. VAN DOMSELAAR.

INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.

RECEIVED: 750601.

BRIEF DESCRIPTION:

    PEIDE ESTIMATES UNKNOWN VARIABLES IN A SYSTEM OF
    FIRST ORDER DIFFERENTIAL EQUATIONS; THE UNKNOWN VARIABLES MAY
    APPEAR NONLINEAR BOTH IN THE DIFFERENTIAL EQUATIONS AND ITS INITIAL
    VALUES; A SET OF OBSERVED VALUES OF SOME COMPONENTS OF THE SOLUTION
    OF THE DIFFERENTIAL EQUATIONS MUST BE GIVEN;

KEYWORDS:

    PARAMETER ESTIMATION,
    DIFFERENTIAL EQUATIONS,
    INITIAL VALUE PROBLEM,
    DATA FITTING.

CALLING SEQUENCE:

    THE HEADING OF THIS PROCEDURE IS:
    "PROCEDURE" PEIDE(N, M, NOBS, NBP, PAR, RV, BP, JTJINV, IN, OUT,
       DERIV, JAC DFDY, JACDFDP, CALL YSTART, DATA, MONITOR);
    "VALUE" N,M,NOBS; "INTEGER" N,M,NOBS,NBP;
    "ARRAY" PAR,RV,JTJINV,IN,OUT; "INTEGER" "ARRAY" BP;
    "PROCEDURE" CALL YSTART,DATA,MONITOR;
    "BOOLEAN" "PROCEDURE" DERIV,JAC DFDY,JAC DFDP;
    "CODE" 34444;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:        <ARITHMETIC EXPRESSION>;
              THE NUMBER OF DIFFERENTIAL EQUATIONS;
    M:        <ARITHMETIC EXPRESSION>;
              THE NUMBER OF UNKNOWN VARIABLES;
    NOBS:     <ARITHMETIC EXPRESSION>;
              THE NUMBER OF OBSERVATIONS; NOBS SHOULD SATISFY  NOBS>=M;
    NBP:      <VARIABLE>;
              ENTRY: THE NUMBER OF BREAK-POINTS; IF NO BREAK-POINTS ARE
                     USED THEN SET NBP=0;
              EXIT:  WITH NORMAL TERMINATION OF THE PROCESS NBP=0;
                     OTHERWISE, IF THE PROCESS HAS BEEN BROKEN OFF (SEE
                     OUT[1]), THE VALUE OF NBP IS THE NUMBER OF BREAK-
                     POINTS USED BEFORE THE PROCESS BROKE OFF;

    PAR:      <ARRAY IDENTIFIER>;
              "ARRAY" PAR[1 : M+NBP];
              ENTRY: PAR[1:M] SHOULD CONTAIN AN INITIAL APPROXIMATION
                     TO THE REQUIRED PARAMETER VECTOR;
              EXIT:  PAR[1:M] CONTAINS THE CALCULATED PARAMETER VECTOR;
                     IF OUT[1]>0 AND NBP>0 THEN PAR[M+1:M+NBP] CONTAINS
                     THE VALUES OF THE NEWLY INTRODUCED PARAMETERS
                     BEFORE THE PROCESS BROKE OFF;
    RV:       <ARRAY IDENTIFIER>;
              "ARRAY" RV[1 : NOBS+NBP];
              EXIT:  RV[1:NOBS] CONTAINS THE RESIDUAL VECTOR AT THE
                     CALCULATED MINIMUM; IF OUT[1]>0 AND NBP>0 THEN
                     RV[NOBS+1:NOBS+NBP] CONTAINS THE ADDITIONAL
                     CONTINUITY REQUIREMENTS AT THE BREAK-POINTS BEFORE
                     THE PROCESS BROKE OFF;
    BP:       <ARRAY IDENTIFIER>;
              "INTEGER" "ARRAY" BP[0 : NBP];
              ENTRY: BP[I], I=1,...,NBP, SHOULD CORRESPOND TO THE INDEX
                     OF THAT TIME OF OBSERVATION WHICH WILL BE USED AS
                     A BREAK-POINT (1<=BP[I]<=NOBS); THE BREAK-POINTS
                     HAVE TO BE ORDERED SUCH THAT BP[I]<=BP[J] IF I<=J;
              EXIT:  WITH NORMAL TERMINATION OF THE PROCESS BP[1:NBP]
                     CONTAINS NO INFORMATION; OTHERWISE, IF OUT[1]>0
                     AND NBP>0 THEN BP[I], I=1,...,NBP, CONTAINS THE
                     INDEX OF THAT TIME OF OBSERVATION WHICH WAS USED
                     AS A BREAK-POINT BEFORE THE PROCESS BROKE OFF;
    JTJINV:   <ARRAY IDENTIFIER>;
              "ARRAY" JTJINV[1 : M, 1 : M];
              EXIT:  THE INVERSE OF THE MATRIX  J' * J  WHERE J DENOTES
                     THE MATRIX OF PARTIAL DERIVATIVES DRV[I] / DPAR[K]
                     (I=1,...,NOBS ; K=1,...,M) AND J' DENOTES THE
                     TRANSPOSE OF J; THIS MATRIX CAN BE USED IF
                     ADDITIONAL INFORMATION ABOUT THE RESULT IS
                     REQUIRED; E.G. STATISTICAL DATA SUCH AS THE
                     COVARIANCE MATRIX, CORRELATION MATRIX AND
                     CONFIDENCE INTERVALS CAN EASILY BE CALCULATED FROM
                     JTJINV AND OUT[2];
    IN:       <ARRAY IDENTIFIER>;
              "ARRAY" IN[0 : 6];
              ENTRY:  IN THIS ARRAY THE USER SHOULD GIVE SOME DATA TO
                      CONTROL THE PROCESS;
              IN[0]:  THE MACHINE PRECISION;
                      FOR THE CYBER 73 A SUITABLE VALUE IS "-14;
              IN[1]:  THE RATIO: THE MINIMAL STEPLENGTH FOR THE
                      INTEGRATION OF THE DIFFERENTIAL EQUATIONS DIVIDED
                      BY THE DISTANCE BETWEEN TWO NEIGHBOURING
                      OBSERVATIONS; MOSTLY, A SUITABLE VALUE IS "-4;
              IN[2]:  THE RELATIVE LOCAL ERROR BOUND FOR THE
                      INTEGRATION PROCESS; THIS VALUE SHOULD SATISFY
                      IN[2]<=IN[3]; THIS PARAMETER CONTROLS THE
                      ACCURACY OF THE NUMERICAL INTEGRATION; MOSTLY,
                      A SUITABLE VALUE IS IN[3]/100;

              IN[3], IN[4]:
                      THE RELATIVE AND THE ABSOLUTE TOLERANCE FOR
                      THE DIFFERENCE BETWEEN THE EUCLIDEAN NORM OF THE
                      ULTIMATE AND PENULTIMATE RESIDUAL VECTOR
                      RESPECTIVELY;
                      THE PROCESS IS TERMINATED IF THE IMPROVEMENT OF
                      THE SUM OF SQUARES IS LESS THAN
                      IN[3] * (SUM OF SQUARES) + IN[4] * IN[4];
                      THESE TOLERANCES SHOULD BE CHOSEN IN ACCORDANCE
                      WITH THE RELATIVE, RESP. ABSOLUTE ERRORS IN THE
                      OBSERVATIONS;
                      NOTE THAT THE EUCLIDEAN NORM OF THE RESIDUAL
                      VECTOR IS DEFINED AS THE SQUARE ROOT OF THE SUM
                      OF SQUARES;
              IN[5]:  THE MAXIMUM NUMBER OF TIMES THAT THE INTEGRATION
                      OF THE DIFFERENTIAL EQUATIONS IS PERFORMED;
              IN[6]:  A STARTING VALUE USED FOR THE RELATION BETWEEN
                      THE GRADIENT AND THE GAUSS-NEWTON DIRECTION (SEE
                      [1]);  IF THE PROBLEM IS WELL CONDITIONED THEN A
                      SUITABLE VALUE FOR IN[6] WILL BE  0.01; IF THE
                      PROBLEM IS ILL CONDITIONED THEN IN[6] SHOULD BE
                      GREATER, BUT THE VALUE OF IN[6] SHOULD SATISFY:
                      IN[0] < IN[6] <= 1/IN[0];
    OUT:      <ARRAY IDENTIFIER>;
              "ARRAY" OUT[1 : 7];
              EXIT :  IN ARRAY  OUT  SOME BY-PRODUCTS ARE DELIVERED;
              OUT[1]: THIS VALUE GIVES INFORMATION ABOUT THE
                      TERMINATION OF THE PROCESS;
                      OUT[1]=0:  NORMAL TERMINATION;
                      IF OUT[1]>0 THEN THE PROCESS HAS BEEN BROKEN OFF
                      AND THIS MAY OCCUR BECAUSE OF THE FOLLOWING
                      REASONS:
                      OUT[1]=1:  THE NUMBER OF INTEGRATIONS PERFORMED
                                 EXCEEDED THE NUMBER GIVEN IN IN[5];
                      OUT[1]=2:  THE DIFFERENTIAL EQUATIONS ARE VERY
                                 NONLINEAR; DURING AN INTEGRATION THE
                                 VALUE OF IN[1] WAS DECREASED BY A
                                 FACTOR 10000 AND IT IS ADVISED TO
                                 DECREASE IN[1], ALTHOUGH THIS WILL
                                 INCREASE COMPUTING TIME;
                      OUT[1]=3:  A CALL OF DERIV DELIVERED THE VALUE
                                 FALSE;
                      OUT[1]=4:  A CALL OF  JAC DFDY  DELIVERED THE
                                 VALUE FALSE;
                      OUT[1]=5:  A CALL OF  JAC DFDP  DELIVERED  THE
                                 VALUE FALSE;
                      OUT[1]=6:  THE PRECISION ASKED FOR CAN NOT BE
                                 ATTAINED; THIS PRECISION IS POSSIBLY
                                 CHOSEN TOO HIGH, RELATIVE TO THE
                                 PRECISION IN WHICH THE RESIDUAL VECTOR
                                 IS CALCULATED (SEE IN[3]);
              OUT[2]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR
                      CALCULATED WITH VALUES OF THE UNKNOWNS DELIVERED;

              OUT[3]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR
                      CALCULATED WITH THE INITIAL VALUES OF THE
                      UNKNOWN VARIABLES;
              OUT[4]: THE NUMBER OF INTEGRATIONS PERFORMED, NEEDED TO
                      OBTAIN THE CALCULATED RESULT; IF OUT[4]=1 AND
                      OUT[1]>0 THEN THE MATRIX JTJINV CAN NOT BE USED;
              OUT[5]: THE MAXIMUM NUMBER OF TIMES THAT THE REQUESTED
                      LOCAL ERROR BOUND WAS EXCEEDED IN ONE
                      INTEGRATION; IF IT IS A LARGE NUMBER, IT MAY BE
                      BETTER TO DECREASE THE VALUE OF IN[1];
              OUT[6]: THE  IMPROVEMENT OF THE EUCLIDEAN NORM OF THE
                      RESIDUAL VECTOR IN THE LAST ITERATION STEP OF THE
                      PROCESS OF MARQUARDT;
              OUT[7]: THE CONDITION NUMBER OF  J' * J , I.E. THE RATIO
                      OF ITS LARGEST TO SMALLEST EIGENVALUES;
    DERIV:    <PROCEDURE IDENTIFIER>;
              THIS PROCEDURE DEFINES THE RIGHT HAND SIDE OF THE
              DIFFERENTIAL EQUATIONS;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "BOOLEAN" "PROCEDURE" DERIV(PAR, Y, T, DF); "VALUE" T;
              "REAL" T; "ARRAY" PAR,Y,DF;
              ENTRY: PAR,Y,T;
                     PAR[1:M] CONTAINS THE CURRENT VALUES OF THE
                     UNKNOWNS AND SHOULD NOT BE ALTERED;
                     Y[1:N] CONTAINS THE SOLUTIONS OF THE DIFFERENTIAL
                     EQUATIONS AT TIME T AND SHOULD NOT BE ALTERED;
              EXIT:  "ARRAY" DF[1 : N];
                     AN ARRAY ELEMENT DF[I] SHOULD CONTAIN THE RIGHT
                     HAND SIDE OF THE I-TH DIFFERENTIAL EQUATION;
              AFTER A SUCCESSFUL CALL OF DERIV, THE BOOLEAN PROCEDURE
              SHOULD DELIVER THE VALUE TRUE;
              HOWEVER, IF DERIV DELIVERS THE VALUE FALSE, THEN THE
              PROCESS IS TERMINATED (SEE OUT[1]);
              HENCE, PROPER PROGRAMMING OF DERIV MAKES IT POSSIBLE TO
              AVOID CALCULATION OF THE RIGHT HAND SIDE WITH VALUES OF
              THE UNKNOWN VARIABLES WHICH CAUSE OVERFLOW IN THE
              COMPUTATION;
    JAC DFDY: <PROCEDURE IDENTIFIER>;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "BOOLEAN" "PROCEDURE" JAC DFDY(PAR, Y, T, FY); "VALUE" T;
              "REAL" T; "ARRAY" PAR,Y,FY;
              ENTRY: PAR,Y,T;
                     SEE DERIV;
              EXIT:  "ARRAY" FY[1 : N,1 : N];
                     AN ARRAY ELEMENT FY[I,J] SHOULD CONTAIN THE
                     PARTIAL DERIVATIVE OF THE RIGHT HAND SIDE OF THE
                     I-TH DIFFERENTIAL EQUATION WITH RESPECT TO Y[J],
                     I.E. DF[I]/DY[J];
              THE BOOLEAN VALUE SHOULD BE ASSIGNED TO THIS PROCEDURE
              IN THE SAME WAY AS IT IS DONE FOR THE VALUE OF DERIV;
    JAC DFDP: <PROCEDURE IDENTIFIER>;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "BOOLEAN" "PROCEDURE" JAC DFDP(PAR, Y, T, FP); "VALUE" T;
              "REAL" T; "ARRAY" PAR,Y,FP;

              ENTRY: PAR,Y,T;
                     SEE DERIV;
              EXIT:  "ARRAY" FP[1 : N,1 : M];
                     AN ARRAY ELEMENT FP[I,J] SHOULD CONTAIN THE
                     PARTIAL DERIVATIVE OF THE RIGHT HAND SIDE OF THE
                     I-TH DIFFERENTIAL EQUATION WITH RESPECT TO PAR[J],
                     I.E. DF[I]/DPAR[J];
              THE BOOLEAN VALUE SHOULD BE ASSIGNED TO THIS PROCEDURE
              IN THE SAME WAY AS IT IS DONE FOR THE VALUE OF DERIV;
    CALL YSTART: <PROCEDURE IDENTIFIER>;
              THIS PROCEDURE DEFINES THE INITIAL VALUES OF THE INITIAL
              VALUE PROBLEM;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "BOOLEAN" "PROCEDURE" CALL YSTART(PAR, Y, YMAX);
              "ARRAY" PAR,Y,YMAX;
              ENTRY: PAR;
                     PAR[1:M] CONTAINS THE CURRENT VALUES OF THE
                     UNKNOWN VARIABLES AND SHOULD NOT BE ALTERED;
              EXIT:  Y,YMAX;
                     Y[1:N] SHOULD CONTAIN THE INITIAL VALUES OF THE
                     CORRESPONDING DIFFERENTIAL EQUATIONS;
                     THE INITIAL VALUES MAY BE FUNCTIONS OF THE UNKNOWN
                     VARIABLES PAR; IN THAT CASE, THE INITIAL VALUES OF
                     DY/DPAR ALSO HAVE TO BE SUPPLIED; NOTE THAT
                     DY[I]/DPAR[J] CORRESPONDS WITH  Y[5*N+J*N+I]
                     (I=1,...,N , J=1,...,M);
                     YMAX[I], I=1,...,N, SHOULD CONTAIN A ROUGH
                     ESTIMATE TO THE MAXIMAL ABSOLUTE VALUE OF Y[I]
                     OVER THE INTEGRATION INTERVAL;
    DATA:     <PROCEDURE IDENTIFIER>;
              THIS PROCEDURE TAKES THE DATA TO FIT INTO THE PROCEDURE
              PEIDE;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "PROCEDURE" DATA(NOBS, TOBS, OBS, COBS); "VALUE" NOBS;
              "INTEGER" NOBS; "ARRAY" TOBS,OBS; "INTEGER" "ARRAY" COBS;
              ENTRY: NOBS;
                     NOBS HAS THE SAME MEANING AS IN PEIDE;
              EXIT:  "ARRAY" TOBS[0 : NOBS];
                     THE ARRAY ELEMENT TOBS[0] SHOULD CONTAIN THE TIME,
                     CORRESPONDING TO THE INITIAL VALUES OF Y GIVEN IN
                     THE PROCEDURE  CALL YSTART; AN ARRAY ELEMENT
                     TOBS[I], 1<=I<=NOBS, SHOULD CONTAIN THE I-TH TIME
                     OF OBSERVATION; THE OBSERVATIONS HAVE TO BE
                     ORDERED SUCH THAT TOBS[I]<=TOBS[J] IF I<=J;
                     "INTEGER" "ARRAY" COBS[1:NOBS];
                     AN ARRAY ELEMENT COBS[I] SHOULD CONTAIN THE
                     COMPONENT OF Y OBSERVED AT TIME TOBS[I]; NOTE THAT
                     1<=COBS[I]<=N;
                     "ARRAY" OBS[1:NOBS];
                     AN ARRAY ELEMENT OBS[I] SHOULD CONTAIN THE
                     OBSERVED VALUE OF THE COMPONENT COBS[I] OF Y AT
                     THE TIME TOBS[I];

    MONITOR:  <PROCEDURE IDENTIFIER>;
              THIS PROCEDURE CAN BE USED TO OBTAIN INFORMATION ABOUT
              THE COURSE OF THE ITERATION PROCESS; IF NO INTERMEDIATE
              RESULTS ARE DESIRED, A DUMMY PROCEDURE SATISFIES;
              THE HEADING OF THIS PROCEDURE SHOULD BE:
              "PROCEDURE" MONITOR(POST,NCOL,NROW,PAR,RV,WEIGHT,NIS);
              "VALUE" POST,NCOL,NROW,WEIGHT,NIS;
              "INTEGER" POST,NCOL,NROW,WEIGHT,NIS; "ARRAY" PAR,RV;
              INSIDE PEIDE, THE PROCEDURE MONITOR IS CALLED AT TWO
              DIFFERENT PLACES AND THIS IS DENOTED BY THE VALUE OF
              POST:
              POST=1:  MONITOR IS CALLED AFTER AN INTEGRATION OF THE
                       DIFFERENTIAL EQUATIONS; AT THIS PLACE ARE
                       AVAILABLE: THE CURRENT VALUES OF THE UNKNOWN
                       VARIABLES PAR[1:NCOL], WHERE NCOL=M+NBP, THE
                       CALCULATED RESIDUAL VECTOR RV[1:NROW], WHERE
                       NROW=NOBS+NBP, AND THE VALUE OF NIS, WHICH IS
                       THE NUMBER OF INTEGRATION STEPS PERFORMED DURING
                       THE SOLUTION OF THE LAST INITIAL VALUE PROBLEM;
              POST=2:  MONITOR IS CALLED BEFORE A MINIMIZATION OF THE
                       EUCLIDEAN NORM OF THE RESIDUAL VECTOR WITH THE
                       PROCEDURE MARQUARDT (SEE SECTION 5.1.3.1.3) IS
                       STARTED; AVAILABLE ARE THE CURRENT VALUES OF
                       PAR[1:NCOL] AND  THE VALUE OF THE WEIGHT, WITH
                       WHICH THE CONTINUITY REQUIREMENTS AT THE BREAK-
                       POINTS ARE ADDED TO THE ORIGINAL LEAST SQUARES
                       PROBLEM.

DATA AND RESULTS: SEE REF[1].

PROCEDURES USED:

    INIVEC = CP31010,
    INIMAT = CP31011,
    MULVEC = CP31020,
    MULROW = CP31021,
    DUPVEC = CP31030,
    DUPMAT = CP31035,
    VECVEC = CP34010,
    MATVEC = CP34011,
    ELMVEC = CP34020,
    SOL = CP34051,
    DEC = CP34300,
    MARQUARDT = CP34440.

REQUIRED CENTRAL MEMORY :
      IN THE BODY OF PEIDE (3 + NBP) * NOBS +
      N * (13 + N + 7 * M + 7 * NBP) ARRAY
      ELEMENTS ARE DECLARED.

METHOD AND PERFORMANCE:
     PEIDE ESTIMATES UNKNOWN VARIABLES IN THE SYSTEM OF DIFFERENTIAL
     EQUATIONS  DY/DT (T, PAR) = F (T, Y, PAR), BY USING A SET OF
     OBSERVED VALUES OF Y; THE UNKNOWN VARIABLES  PAR  ARE OBTAINED IN
     THE LEAST SQUARES SENSE; AN ELEMENT OF THE RESIDUAL VECTOR IS
     DEFINED BY THE CALCULATED VALUE OF Y MINUS ITS OBSERVED VALUE;
     THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR IS MINIMIZED BY THE
     ITERATION PROCESS OF MARQUARDT; THE DIFFERENTIAL EQUATIONS ARE
     SOLVED BY THE INTEGRATION PROCESS OF GEAR; A MULTIPLE SHOOTING
     TECHNIQUE HAS BEEN IMPLEMENTED TO IMPROVE BAD STARTING VALUES OF
     THE UNKNOWNS; IF THIS TECHNIQUE IS USED, ONE HAS TO GIVE SOME
     BREAK-POINTS, I.E. TIMES OF OBSERVATIONS WHERE A NEW INITIAL
     VALUE PROBLEM SHOULD BE STARTED; THE NEW INITIAL VALUES OF Y
     BECOME EXTRA UNKNOWN VARIABLES AND THE CONTINUITY REQUIREMENTS
     AT THE BREAK-POINTS ARE ADDED WITH SOME WEIGHTING FACTOR TO THE
     LEAST SQUARES PROBLEM; FOR DETAILS SEE REF[1].

REFERENCES:
     [1]: B. VAN DOMSELAAR,
          NONLINEAR PARAMETER ESTIMATION IN INITIAL VALUE PROBLEMS,
          MATH. CENTRE, AMSTERDAM (TO APPEAR).

EXAMPLE OF USE:

     THE PARAMETERS PAR[1:3] IN THE DIFFERENTIAL EQUATIONS
        DY[1]/DT = - (1 - Y[2]) * Y[1] + EXP(PAR[2]) * Y[2],
        DY[2]/DT = EXP(PAR[1]) * ((1 - Y[2]) * Y[1] - (EXP(PAR[2])+
                   +EXP(PAR[3])) * Y[2]),
     WITH 23 OBSERVATIONS OF Y[2], MAY BE OBTAINED BY THE FOLLOWING
     PROGRAM, THAT CONSISTS OF
          1: A CODE PROCEDURE WHICH TAKES CARE OF THE OUTPUT OF THE
          EXAMPLE PROGRAM. IT ALSO INTERPRETS THE NUMERICAL DATA
          THAT CAN BE USED TO OBTAIN STATISTICAL RESULTS;
          2: THE USERS PROGRAM IN WHICH THE PROBLEM EXAMPLE IS DEFINED.

"CODE" 34445;

THE USER PROGRAM READS:

"BEGIN" "INTEGER" I,M,N,NOBS,NBP; "REAL" TIME,FA;
 "ARRAY" PAR[1:6],RES[1:26],JTJINV[1:3,1:3],IN[0:6],OUT[1:7];
 "INTEGER" "ARRAY" BP[0:3];

 "PROCEDURE" COMMUNICATION(POST,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV,
                           IN,OUT,WEIGHT,NIS);
 "VALUE" POST,FA,N,M,NOBS,NBP,WEIGHT,NIS;
 "INTEGER" POST,N,M,NOBS,NBP,WEIGHT,NIS; "REAL" FA;
 "ARRAY" PAR,RES,JTJINV,IN,OUT; "INTEGER""ARRAY" BP;
 "CODE" 34445;

 "BOOLEAN" "PROCEDURE" JAC DFDP(PAR,Y,X,FP);
    "REAL" X; "ARRAY" PAR,Y,FP;
    "BEGIN" "REAL" Y2; Y2:=Y[2];
       FP[1,1]:=FP[1,3]:=0;
       FP[1,2]:=Y2*EXP(PAR[2]);
       FP[2,1]:=EXP(PAR[1])*(Y[1]*(1-Y2)-(EXP(PAR[2])+EXP(PAR[3]))*Y2);
       FP[2,2]:=-EXP(PAR[1]+PAR[2])*Y2;
       FP[2,3]:=-EXP(PAR[1]+PAR[3])*Y2;
       JAC DFDP:="TRUE"
    "END" JAC DFDP

 "PROCEDURE" DATA(NOBS,TOBS,OBS,COBS);
    "VALUE" NOBS; "INTEGER" NOBS;
    "ARRAY" TOBS,OBS; "INTEGER" "ARRAY" COBS;
    "BEGIN" "INTEGER" I;
       TOBS[0]:=0;
       OUTPUT(61,"("*,4/,4B,"("THE OBSERVATIONS WERE:")",
       //,B,"("I")",3B,"("TOBS[I]")",3B,"("COBS[I]")",3B,
       "("OBS[I]")",/")");
       "FOR" I:=1 "STEP" 1 "UNTIL" NOBS "DO"
       "BEGIN" INREAL(70, TOBS[I]); ININTEGER(70, COBS[I]);
          INREAL(70, OBS[I]);
          OUTPUT(61,"("ZD,3B,ZD.4D,6B,D,6B,.4D,/")",I,TOBS[I],COBS[I],
          OBS[I])
       "END"
    "END" DATA;

 "PROCEDURE" CALL YSTART(PAR,Y,YMAX);
    "ARRAY" PAR,Y,YMAX;
    "BEGIN" Y[1]:=YMAX[1]:=YMAX[2]:=1;
       Y[2]:=0
    "END" CALL YSTART;

 "BOOLEAN" "PROCEDURE" DERIV(PAR,Y,X,DF);
    "REAL" X; "ARRAY" PAR,Y,DF;
    "BEGIN" "REAL" Y2; Y2:=Y[2];
       DF[1]:=-(1-Y2)*Y[1]+EXP(PAR[2])*Y2;
       DF[2]:=EXP(PAR[1])*((1-Y2)*Y[1]-(EXP(PAR[2])+EXP(PAR[3]))*Y2);
       DERIV:="TRUE"
    "END" DERIV;

 "BOOLEAN" "PROCEDURE" JAC DFDY(PAR,Y,X,FY);
    "REAL" X; "ARRAY" PAR,Y,FY;
    "BEGIN" FY[1,1]:=-1+Y[2];
       FY[1,2]:=EXP(PAR[2])+Y[1];
       FY[2,1]:=EXP(PAR[1])*(1-Y[2]);
       FY[2,2]:=-EXP(PAR[1])*(EXP(PAR[2])+EXP(PAR[3])+Y[1]);
       JAC DFDY:="TRUE"
    "END" JAC DFDY;

 "PROCEDURE" MONITOR(POST,NCOL,NROW,PAR,RES,WEIGHT,NIS);
  "VALUE" POST,NCOL,NROW,WEIGHT,NIS;
  "INTEGER" POST,NCOL,NROW,WEIGHT,NIS; "ARRAY" PAR,RES;;

 OUTPUT(61,"("2/,30B,"("E S C E P - PROBLEM")",3/")");
 M:= 3; N:=2; NOBS:=23; NBP:=3;
 PAR[1]:=LN(1600); PAR[2]:=LN(.8); PAR[3]:=LN(1.2); IN[0]:="-14;
 IN[3]:="-4; IN[4]:="-4; IN[5]:=50; IN[6]:="-2;
 IN[1]:="-4; IN[2]:="-5;
 BP[1]:=17; BP[2]:=19; BP[3]:=21;
 FA:=4.94;
 "COMMENT"  FA  DENOTES THE ALPHA-POINT OF THE FISHER-DISTRIBUTION;

 COMMUNICATION(1,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,0,0);
 TIME:=CLOCK;

 PEIDE(N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,DERIV,JAC DFDY,JAC DFDP,
 CALL YSTART,DATA,MONITOR);

 TIME:=CLOCK-TIME;
 COMMUNICATION(2,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,0,0);
 OUTPUT(61,"("3/,5B,
 "("THE CALCULATION IN PEIDE CONSUMED")",B,ZZD.DD,2B,
 "("SECONDS")",*")",TIME)
"END"

 THIS PROGRAM DELIVERS:

                              E S C E P - PROBLEM

          STARTING VALUES OF THE PARAMETERS
                    +.7377759" +1
                    -.2231436" +0
                    +.1823216" +0

NUMBER OF EQUATIONS   : 2
NUMBER OF OBSERVATIONS:23

MACHINE PRECISION                              :+.1"-13
RELATIVE LOCAL ERROR BOUND FOR INTEGRATION     :+.1" -4
RELATIVE TOLERANCE FOR RESIDUE                 :+.10" -3
ABSOLUTE TOLERANCE FOR RESIDUE                 :+.10" -3
MAXIMUM NUMBER OF INTEGRATIONS TO PERFORM      : 50
RELATIVE STARTING VALUE OF LAMBDA              :+.10" -1
RELATIVE MINIMAL STEPLENGTH                    :+.10" -3

BREAK-POINTS ARE THE OBSERVATIONS : 17  19  21

THE ALPHA-POINT OF THE F-DISTIBUTION : 4.94

    THE OBSERVATIONS WERE:

 I   TOBS[I]   COBS[I]   OBS[I]
 1    0.0002      2      .1648
 2    0.0004      2      .2753
 3    0.0006      2      .3493
 4    0.0008      2      .3990
 5    0.0010      2      .4322
 6    0.0012      2      .4545
 7    0.0014      2      .4695
 8    0.0016      2      .4795
 9    0.0018      2      .4862
10    0.0020      2      .4907
11    0.0200      2      .4999
12    0.0400      2      .4998
13    0.0600      2      .4998
14    0.0800      2      .4998
15    0.1000      2      .4998
16    1.0000      2      .4986
17    2.0000      2      .4973
18    5.0000      2      .4936
19   10.0000      2      .4872
20   15.0000      2      .4808
21   20.0000      2      .4743
22   25.0000      2      .4677
23   30.0000      2      .4610

NORMAL TERMINATION OF THE PROCESS

LAST INTEGRATION WAS PERFORMED WITHOUT BREAK-POINTS

EUCL. NORM OF THE LAST RESIDUAL VECTOR :.1430776" -3
EUCL. NORM OF THE FIRST RESIDUAL VECTOR:.1331071" +1
NUMBER OF INTEGRATIONS PERFORMED       : 12
LAST IMPROVEMENT OF THE EUCLIDEAN NORM :.2223694" -4
CONDITON NUMBER OF  J'*J               :.2582882" +3
LOCAL ERROR BOUND WAS EXCEEDED (MAXIM.): 37

 PARAMETERS            CONFIDENCE INTERVAL
+.6907670" +1            +.3209313" -3
-.1003941" -1            +.1687774" -3
-.4605292" +1            +.1942501" -2

CORRELATION MATRIX           COVARIANCE MATRIX
                             +.6949857" -8 +.1407628" -8 -.9129848" -8
+.3851320" +0                              +.1922119" -8 -.1414245" -7
-.2170393" +0 -.6392889" +0                              +.2546094" -6

          THE LAST RESIDUAL VECTOR

               I    RES[I]
               1  +.1748" -5
               2  -.2905" -4
               3  +.2814" -4
               4  -.3879" -4
               5  +.3069" -4
               6  +.3101" -4
               7  -.2019" -4
               8  -.3887" -5
               9  +.1052" -4
              10  +.1391" -4
              11  -.5109" -4
              12  +.2384" -4
              13  -.1156" -5
              14  -.2616" -4
              15  -.5116" -4
              16  +.2244" -4
              17  +.6794" -4
              18  -.1418" -4
              19  +.2087" -4
              20  -.1980" -4
              21  -.3476" -4
              22  -.2245" -4
              23  +.1886" -4

     THE CALCULATION IN PEIDE CONSUMED 108.57  SECONDS

SOURCE TEXT(S):
"CODE" 34444;