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;