- Email: [email protected]
Program to estimate parameters of linear systems without numerical differentiation
Computer Programs in Biomedicine 8 (1978) 135-140 © Elsevier/North-Holland Biomedical Press
P R O G R A M TO E S T I M A T E P A R A M E T E R S OF L I N E A R SYSTEMS WITHOUT N U M E R I C A L DIFFERENTIATION B. KANY~,R and J. ERODI Semmelweis University of Medicine, Computing Group, 1089 Budapest, Hungary This paper describes a computer program for estimating the parameters of a linear differential equation system with constant coefficients by use of a nonlinear least-squares method. For minimization the sum of squares of an existing standard program, the Gauss-Newton gradient procedure, is employed. The differential equation system is solved by the Taylor expansion method. The advantage of this approach is that the derivatives with respect to the parameters are available without numerical differentiation. Therefore the inaccuracy inherent in numerical differentiation and the problem of choosing the modification of the parameters are eliminated. The given procedure is applicable for all the first order gradient methods. The presented method was tested with generated data from a four-compartmental model. Nonlinear parameter estimation
Gradientmethod
Compartment
1. Introduction
Differentialequations
binatior~ of a Runge-Kutta and the Adams-Moulton scheme. In case of a linear system with constant coefficients an eigenvalue-eigenvector method can also be used. The procedure given by Jennrich and Bright [6] expresses the derivatives o f y with respect to the parameters akl in terms of the eigenvectors and eigenvalues. In this way they avoid numerical derivating and the computational errors stem only from the solution of the eigenvalue problem. This paper described a procedure, which uses a gradient method to fit a linear system with constant coefficients, but to avoid the numerical differentiation the partials are formed by using the Taylor-expansion method. A numerical approximation is only used at the finite Taylor expansion. The procedure is applicable for all the gradient methods which need the first order partials (e.g., steepest descent, Gauss-Newton and Marquardt methods [1]).
In tracer kinetic investigations the systems studied can often be described by linear systems of differential equations with constant coefficients. In vector matrix form they can be written: dy -d-[ = A y (1) According to the method of least squares, the optimum value of the parameters a u may be obtained by minimizing the sum of squares of deviations between the experimental and the calculated data. To minimize the sum of squares, direct-search and gradient-search methods [ 1,8] are applicable, When the gradient methods are applied the greatest problem is to form the partial derivatives ~y/aakl. Berman et al. [2] calculate the partial derivatives with respect to the parameters akt numerically, changing the parameters by 1-2%. This is a generally applicable method but the approximation of the partials may be less accurate or the exact calculation requires much computer time. In the program by Buell et al. [3] the partials are formed during the solution of the system, To minimize the sum of squares they used the NewtonRaphson gradient method. This method requires also the second order partial derivatives. The program integrates the differential equation system by acom-
2. Mathematical background As mentioned above the greatest problem of fitting the linear model (1) and to minimize the sum of squares by gradient methods is to form the ay/~akl partial derivatives, where y is the state variable vector and a u are the elements of matrix A, and these are the parameters. 135
B. Kanytir, J. ErOdL Parameters of linear systems without numerical differentiation
136
The vector y and the ~y/aakt vectors as well, can be achieved by using stepwise the Taylor expansion with the stepsize h: N
y(t + h) ~ ~ m=o
(hA)m m! "y(t)
(2)
and by partial differentiation
ay(t +h) ~ ing=0 (hA)mq ay(t) I N~ h m aA m .ly aakt m! ] ~ + m=l-~.. ~ (t) -
(3) (It can be proven that at N = 4 the Taylor-expansion is equivalent to the fourth-order Runge-Kutta method.) According to the rules of matrix analysis we can give a rather convenient form [7] to calculate the derivatives of the matrix Am:
aA m ~akt Ibt~-l) ' b~r~-r) btrg-1) " bt~-r) ... btrk-t)" bt~ -r) b~rk_l) " btr~_r) b~rk_l) " btr~_r).., b~rk_l) " btsm_r) m
estimated approximately by the formula: ( N - 1)s 3 + K(p + 1)(N+ 2)s 2 where s is the dimension of matrix A and p the number of the parameters. The accuracy of the presented procedure to calculate the partial derivatives ay/~akt was investigated by use of the model given in the sample run. For comparison of the different methods this study has investigated both the one-sided and the two-sided numerical differentiations. The results of this test indicate that the one-sided approximation may differ several times from the correct values. The partial derivatives got from the two-sided numerical differentiation and from the presented method agreed with the correct values well (within 5%). But the numerical differentiation was mostly less accurate than our method. The onesided numerical differentiation needs about the same, the two-sided about twice as many operations as the presented procedure. Even if comparisons between the computer time needed are difficult - if they do not apply to the very same problem on the same comp u t e r - i t couldbe mentionedthat, e.g.,the program given by Buell et al. [3] uses about three times more CPU time than our program.
o,.
r=l Lb~ -1) " bt'~ -r)
b(srk-l) " bt~ -r) ... b(~-l) " btm-r)
3. Program description
(4) where b ~ - l ) is the element of the ith row and kth column of the matrix A r-1. The symbol s represents the number of equations in the system, the dimension of matrix A.
To minimize the sum of squares the program uses (with suitable modification according to our problem) the Gauss-Newton method given by Jennrich and Sampson in program BMDX85 [5]. Figure 1 gives the
For testing the accuracy of the Taylor-series method the error o f y may be approximated by:
where K is the number of steps in the integration of the system [7]. The advantage of the expression (5) is that only few calculations are required, but according to our experience it is a severe condition, therefore if it fails we use the method of halving the stepsize. The number of multiplications necessary to solve the system and to compute the ay/aakt partials can be
1
~'~[
1 STEP
[ I
PLOT
] I
Fig. 1. Structure of the program.
137
B. Kany~r, Z Erfdi, Parame~ of linear sys~ms without num~ical diffemnt~t~n PROBLEM IDENTIF.
ROSE-BENg./TIST
NuM§ER OF THE PARAMETERS
6
NUMBER OF THE OBSERVATIONs
19
NUMOER OF THE iNDEP,VAR|ABLES
1
NUMBER OF THE ~EPENCE NT VARIAB,
2
HuMB|R OF THE dEIGHT
|
VARIAB,
VE|GHT CODE
"1 *1
NUMBER OF THE O I F F . E ~ * °S
4
NUMBER OF TkE ?RIRER PAR,oS
4
1NP~T OATS FORMAT
|
I
(ISFS,I)
NUMBER OF MAXIMUM FUNCTION EVAL,
21
ToLEB,OF FUNCTIDM CHANG|NG
D,|OIEBe
T~LER,CR|TER, AT ~ATR*INvERS,
Q,BIII|I
N~MBER OF THE TAYLOR POWER
4
TA START,VALUE OF IkOEP,VAR
Q,P
ERRCR LIMIT AT soL,OF DoE,S,
Bl*illl
STARTING VALUES oF THE STATE VAR.S INOEcES OF NON-ZERO NATR,ELES
I
lO|,l I 2
2 1
I*l ) 2
O,O
0.0
4 3
F~.2a. Problemdescription on output.
structure of the program with the following subroutines: MAIN The purpose of MAIN is to read input data and print the description of the problem (fig. 2a). Input data Number of the primary and secondary parameters, number of the differential equations, number of data points, number of curves to be fitted simultaneously, number of weight factors, weight codes in case of weighting by measured data, tolerance at nonlinear regression, initial values of the system, error limit in solution of the system, structure of the matrix A given by the indices of the nonzero matrix elements, data format, data, limits and guesses of the parameter values,
MATCO The MATCO subroutine computes the powers of matrix and the aAm/Oakl partials given in eq. (4).
SETDIF It solves the differential equation system and calculates the ay/aakt partial derivatives according to the formulas (2) and (3). Testing of the accuracy is made by relation (5) and by halving of stepsize. OBSEQ The observation equations are defined by supplied FORTRAN statements. The routine OBSEQ describes the connection of the model and response function, e.g., measuring the sum of two compartments.
B. Kany6r, .1..Er6di, Parameters o/linear systems without numerical differentiation
138
MINIMA
lEO
MAXIM~
l,idO|£
IT£RITIo~
ERROR M£AN S@uARE
0 O
7.6138E
O.g O0
l.Og|~E
0.0
O,O
L.OOegE OJ
l.Ji|OE
O]
03
3.d|OdE
03
5.~llE
ld
6.|RIOE 5.8839E
It El
1.0~0E-03 7,4635E-04
2.|~00E-|$ 2.3283E°10
2.0704E012 4.5612E°12
1.1462E002 2.)621E°E2
t,3543E-12 3.7336E-12
203283E-10 2,3283E-10
4.9904E | l 4,2928E O~ 4.5194[ Id
1.3275E-Et 6.6374EoR2 2.9~43E-02
6.0572E-02 6.9665E-02 l.geEtE-Ot
1.8668E002 l.i2d4E-|2 5,2714E-13
L,ES96E-02 1.7dlgE-|2 4,|2| E-E2
3.9SOLE 3,44|1E 2.O~64E
7.8SO6t-i2 ?.9567Eol2
2.1709E Oi
9.8740E OR
2,ill6E 2,OIllE
tO eO
9.9913E 9.9~05E
II O~
00 OI 0~
9.99|6E 9.99|6E 9.9906E 9.9916E 9.9906E 9.9906E
|d Ol J~ go El OI
l.S4S9E 6.3447E
7
O
2.3467E i |
2,4295E-E2
l.lO?9E-IL
7.|OIEE-O3
8 9
I |
3,4838E'13 9.8768E'14
2.5017E-12 2.5096E-O2
I.IIlIE-II l.ig16E-01
7.9699E-03 8*0262E-R3
II 1! 12 13 14
I ] 2
9.6507E'04 9.8506E*E4 9.8$1LE'O4
4
9.8529E°|4
2 4
9,8506E'04 9.8539E°04
2.5~91E-02 2.5096E-E2 2.5096E-|2 2.5g96E-R2 2.5~96E°02 2.5~96E-12
l,iil6E-II I.iII6E-iI I.OgI6E*I& I.I|L6E°O! t.ROI6EoOt 1.0|16£-0!
1S
02
5.illiE-g3 ~.1133E0|3
1 I
4 5 6
5.~OIOE
l.leOiE-02 1.|~85Eo02
| 1 1
2 )
OQ
0.|
PARAMETERS
7.5269E 03 I . t ~ 7 9 E 13 4 . 0 4 L 4 [ 02 3 . $ 6 6 9 E 02
1
O.O
12 It
7.8184E-R2
il 00 iO
6.R264E-03
7.6217E-R2
8.02~4E-03
?.0217E-R3
8.|2~4E-|3
7.8217E*02 ?.8217E-E2 7.8217E°|2 7.6217E-E2
2.1NllE 2.0Site 2.o~ltE 2.OOlIE 2.OItlE 2.lillE
T,4124Eo|4
1.6867E-03
6.0264E0|3 8.|264E°03 8.|264E-|3
oo IO OI
1.5259E005 6.2842E O~ 9.eII8E 9.7087E I.|151E
id O| O~
ASYMPTOTIC S~A~OARO OEVIATIONS OF THE PARAMETERS 4.9351EoE5
9.8151E-05
i,~483E-OS
4.41eSE-SJ
ASYMPTOTIC CORRELATIo N MATRIX oF THE PARAHETERS I
t 2
2
3
4
5
6
t.OO~il
0,79672
0*66579
00.77737
.0039473
°0.50638
3
0,79672 1.66879
l,liiJi e*417|1
0,41711 l,gliiO
od.6997g -~.667~4
-0,48t19 0E,06797
e0.72202 0~,43154
4 S
°0.77737 -i,39473
°|,067|4 °0.16797
0,10827
0.78647 *i,07?73
6
-i.SI638
*0.69970 -i*4Oll9 "1*72212
-i043154
L.Iilig 1018827 #078647
|,lifll -I,i7773
I,lilll
ESTIMATED MATRIX OF THE OIFF,EQeSYSTEM
o 8 . 1 1 J t 6 l • ! 2 S 1| d,|0~16 -|,0531~
~,i
•e
CASE
I. 0 I,O
l" ! O,O
e . e 78~ 2
o. 0
0. IIdSg3 °11,07822 e. |
RESIDUALS
Fig. 2b. The error me a n square and t he pa ra ma t e r values after
I,E
each iteration, t he standard deviations, correlations and the estimated parameter values after the last iteration.
iNPUT VAR|AGLES
53LUT|oN OF THE [IFF.EQ.SYSTEH
1 2 ] 4 5
i.|46 4.229 g.$62 8.3g9 0,49L
-i,176 -E.294 -|,528 5.223 0,043
l,ill 2.J16 ).14| 4,006 5.|00
916.111 623.481 7S|,glg 605.001 6200000
2752.w4| 3416.1Ja 4Beo,eoI 4520.cd8 4067.~d|
1 2 3 4 ~
l,l~l 2.Jd9 3.|8| 4,081 S,088
6 ? 6 9 18 tl 12 13 14
0.112 -|.|$2 -0.002 -8,260 -0.363 -d.672 -8.541 -8.411 0.273
-0,906 *i,227 "4.145 2.261 - e.688 2,?73 1,160 -2,|62 1,320
8.|06 10.000 15.000 28,J0| 380065 40.000 50.J8| 6~,|06 78,g00
493.480 42?.g80 318.000 257,006
6 7 e 9 le 11 12 12 14
8.880 10.|80 l$.Odo 20.odo
178,|0J 163.000 152.800 1430880
6016,~80 6StA.~O| 7276.~86 7614.e | 8 7644,~08 73660088 6986.~80 6584.~ 0| 6|96.~88
15 16 17 18 19
8.132 8.376 8.095 0.333 0.123
°1.778 -4.J27 4.270 L.633 -~,038
6JeLJO 900088 180*|0| 118.860 12~0696
134°100 126.088 118.000 111.888 194,060
5810.~86 5461.~8| 3135.~00 4821.~08 4316,~8
15 16 17 18 19
202.008
F ~ . 2e. Residuals and input data.
9~,588 U2.277 74.944 68.469 62.751 49,209
9.~76 17.576 24,736 3Q,969 36.436 48.860
1.93T i.z4e 8.294 0.488 0.712 1,493
J.~11 i.,j8 8,42d 0,~53 f,181 0.358
48.080 5J,JiJ 68.080 78,080
42.785 31.000 25.727 20.236 17.087 16.354 1~.241 14.273
54.611 63.129 66.587 66.676 63.852 60,319 560734 53.260
2.049 3.372 4.450 $~789 6.263 6.320 6.136 5.656
0.635 1.699 3.237 T.~96 12,~57 17.~07 21,~86 26.502
~.OdO 90.000 18~.88J 110.086 12~.868
13.387 12.562 L1,790 11,867 16.367
58.029 460962 44.062 41.377 36,036
5.542 50223 4.913 4.616 4,335
31.~41 35.251 39.215 42.940 46.489
3~.88!
F ~ . 2d. S ol ut i on of t he differential e q u a t i o n system a ~ e r the Mst iteration.
139
B. Kanydr, J. Erbdi, Parameters of linear systems without numerical differentiation PLOTS OF TME tEIIOUAL$ I*0109 2*Oil| 3*6J19 4°6166 5*1|9d 0o0990 19*|910 IS,OO|9 2 29*0g90
IO|$EmV,oPtIOICT.wg.) a
| a
)9*6099
114 YM! INO|P|UDtNT VaRIAOk| i | j | t ~ 1 I t I I 2
I
l
4|*1||9
1
5906990
2
I
60.6010
|
2
t
70,9600
I
6900999
~
90.0018
2
I
2
t
1J|,||90
I
119.6990
laO*60dO
|
2
t
2 I -5.d352
l *4.0094
1 -2.96)6
I
I
"1.9570
6 PUg~kF~|~ 4EO|ZTAS AtJ~VA I v l Z S I I N T E S N | I
*0.9320
2
t I
I
0.0937
1.1195
I 2.1453
I
t
).1711
4,1969
I S.22~7
0007d4 : t
Fig. 29. Plots of the residuaL. MINIZ and STEP These subroutines minimize the sum of squares, estimates the parameters, the standard deviations and the correlation matrix by use of the GaussNewton method. They are only slightly modified from the BMDX85 subroutines. MINIZ prints the values of the parameters and the solution of the differ°ntial equation system.
According to the output code the predicted values, the observed data and the residuals may be printed and plotted. The subroutine PLOT from IBM SSP library plots the defined data.
4. Sample run The compartment model used for this run is shown in fig. 3. The model is used to analyse the rose-bengal kinetics in the investigation of liver function [9]. That kind of model is used in our experimental investigations too. The tracer is injected into compartment 1. The rate equations of the model in matrix form are the following:
dy _i Or21
"--(O~12+ (x32)
0
cza2 0
--~43
~ -
•y [0
~x4a
B. Kanydr, J. Er6di, Parameters of linear systems without numerical differentiation
140 ( ~ - - ~ ( ~
"3___~2~
~,___.~3~
~,=I0 ~,,~~23 ~
/2~
f2
Fig. 3. Compartmental model used for testing the program. The three angles represent observations that are directly proportional to the contents of the compartments, the proportional coefficients are o t 1, o21, o22 and a 23. In case of our rose-bengal investigations two data sets fl(t i) and f2(ti) are to be fitted simultaneously. where y is the amount of tracer in the different compartments and akt are the transport coefficients, the primary parameters. The three angles represent observations and the secondary parameter oij is the fraction of compartment ] measured by detector i. To test the program simulated data sets (fl(ti) and f2(ti) ) were generated. The exact values of the data were rounded to three decimal digits. The values of the primary parameters were: al2 = 0.025
Or21 = 0.100
(~32 = 0.008
0t43 = 0.080
and the secondary parameters, the proportional coefficients o21 = 2 and 022 = 023 = 10. The proportional coefficient o1 z was considered as a constant value and not as a parameter, ozl = 10. Figures 2 a - e show the result of the run. The problem and model description is given in fig. 2a. The weight code (q) determines the weighting factors (wi) in the sum of squares, namely wi = f(ti) q. The next figure (2b) shows the error mean square and the parameter values after each iteration, the standard deviations, correlations and the estimated parameters in matrix form after the last iteration. The residuals and the input data, the variables time, fl(ti) and f2 (ti) are printed in fig. 2c. The program gives the solution of the differential equations (fig. 2d) and the plot of the residuals (fig. 2e), after the last iteration, According to the last iteration the estimated values of the parameters approach the original values quite well.
5. Hardware and software specifications
was tested on an IBM 360/75 system. In the case of a three-compartment model with four parameters [4] the total t processor h e gtime parameters. u wase 5 - 1s5 ssdepending e s on initial of the The new version of the program runs on an RIAD20 with DOS. Input unit is the card reader and the output is given at the printer. The required memory is 100 kbyte in case of 12 differential equations 15 parameters and 7 value of the power of Taylor-serie. DOS requires 20 kbyte. The sample run with 4 equations, 6 parameters and two data sets uses 10 min total time (read input, processing and print output).
6. Mode of availability The FORTRAN IV program list may be obtained from the authors at the following address: Semmelweis University of Medicine, Computer Group, 1089 Budapest, Kulich Gy. 5, Hungary.
References [1 ] Y. Bard, Nonlinear Parameter Estimation (Academic Press, New York, London, 1974). [2] M. Berman and M.F. Weiss, SAAM Manual (National Institutes of Health, Bethesda, 1971). [3] J. Buell, R. Kalaba, A. Yakush and E. Ruspini, Comput.
Prog. Biomed. 2 (1971) 8. [4] S. Gy~Srgyiand B. Kany~r, Acta Biochem. Biophys. (Acad. Sci. Hung) 7 (1972) 359. [5 ] R.I. Jennricb and P. Sampson, in: Biomedical Computer Programs (X-Series suppl.) ed. W.J. Dixon (Univ. Calif.
Press, 1972). [6] R.I. Jennrich and P.B. Bright, Technometrics 17 (1975) 245. [7] B. Kany~irand J. T6th, Alkalmazott Matematikai Lapok (Letters of Applied Mathematics) in Hungarian (1978) in press. [8] w. Murray, Numerical Methods for Unconstrained Optimization (Academic Press, New York, London, 1972). 19] A.D. Waxman, P.A. Leins and J.K. Siemsen, Comput. Biomed. Res. 5 (1973) 1.
P R O G R A M TO E S T I M A T E P A R A M E T E R S OF L I N E A R SYSTEMS WITHOUT N U M E R I C A L DIFFERENTIATION B. KANY~,R and J. ERODI Semmelweis University of Medicine, Computing Group, 1089 Budapest, Hungary This paper describes a computer program for estimating the parameters of a linear differential equation system with constant coefficients by use of a nonlinear least-squares method. For minimization the sum of squares of an existing standard program, the Gauss-Newton gradient procedure, is employed. The differential equation system is solved by the Taylor expansion method. The advantage of this approach is that the derivatives with respect to the parameters are available without numerical differentiation. Therefore the inaccuracy inherent in numerical differentiation and the problem of choosing the modification of the parameters are eliminated. The given procedure is applicable for all the first order gradient methods. The presented method was tested with generated data from a four-compartmental model. Nonlinear parameter estimation
Gradientmethod
Compartment
1. Introduction
Differentialequations
binatior~ of a Runge-Kutta and the Adams-Moulton scheme. In case of a linear system with constant coefficients an eigenvalue-eigenvector method can also be used. The procedure given by Jennrich and Bright [6] expresses the derivatives o f y with respect to the parameters akl in terms of the eigenvectors and eigenvalues. In this way they avoid numerical derivating and the computational errors stem only from the solution of the eigenvalue problem. This paper described a procedure, which uses a gradient method to fit a linear system with constant coefficients, but to avoid the numerical differentiation the partials are formed by using the Taylor-expansion method. A numerical approximation is only used at the finite Taylor expansion. The procedure is applicable for all the gradient methods which need the first order partials (e.g., steepest descent, Gauss-Newton and Marquardt methods [1]).
In tracer kinetic investigations the systems studied can often be described by linear systems of differential equations with constant coefficients. In vector matrix form they can be written: dy -d-[ = A y (1) According to the method of least squares, the optimum value of the parameters a u may be obtained by minimizing the sum of squares of deviations between the experimental and the calculated data. To minimize the sum of squares, direct-search and gradient-search methods [ 1,8] are applicable, When the gradient methods are applied the greatest problem is to form the partial derivatives ~y/aakl. Berman et al. [2] calculate the partial derivatives with respect to the parameters akt numerically, changing the parameters by 1-2%. This is a generally applicable method but the approximation of the partials may be less accurate or the exact calculation requires much computer time. In the program by Buell et al. [3] the partials are formed during the solution of the system, To minimize the sum of squares they used the NewtonRaphson gradient method. This method requires also the second order partial derivatives. The program integrates the differential equation system by acom-
2. Mathematical background As mentioned above the greatest problem of fitting the linear model (1) and to minimize the sum of squares by gradient methods is to form the ay/~akl partial derivatives, where y is the state variable vector and a u are the elements of matrix A, and these are the parameters. 135
B. Kanytir, J. ErOdL Parameters of linear systems without numerical differentiation
136
The vector y and the ~y/aakt vectors as well, can be achieved by using stepwise the Taylor expansion with the stepsize h: N
y(t + h) ~ ~ m=o
(hA)m m! "y(t)
(2)
and by partial differentiation
ay(t +h) ~ ing=0 (hA)mq ay(t) I N~ h m aA m .ly aakt m! ] ~ + m=l-~.. ~ (t) -
(3) (It can be proven that at N = 4 the Taylor-expansion is equivalent to the fourth-order Runge-Kutta method.) According to the rules of matrix analysis we can give a rather convenient form [7] to calculate the derivatives of the matrix Am:
aA m ~akt Ibt~-l) ' b~r~-r) btrg-1) " bt~-r) ... btrk-t)" bt~ -r) b~rk_l) " btr~_r) b~rk_l) " btr~_r).., b~rk_l) " btsm_r) m
estimated approximately by the formula: ( N - 1)s 3 + K(p + 1)(N+ 2)s 2 where s is the dimension of matrix A and p the number of the parameters. The accuracy of the presented procedure to calculate the partial derivatives ay/~akt was investigated by use of the model given in the sample run. For comparison of the different methods this study has investigated both the one-sided and the two-sided numerical differentiations. The results of this test indicate that the one-sided approximation may differ several times from the correct values. The partial derivatives got from the two-sided numerical differentiation and from the presented method agreed with the correct values well (within 5%). But the numerical differentiation was mostly less accurate than our method. The onesided numerical differentiation needs about the same, the two-sided about twice as many operations as the presented procedure. Even if comparisons between the computer time needed are difficult - if they do not apply to the very same problem on the same comp u t e r - i t couldbe mentionedthat, e.g.,the program given by Buell et al. [3] uses about three times more CPU time than our program.
o,.
r=l Lb~ -1) " bt'~ -r)
b(srk-l) " bt~ -r) ... b(~-l) " btm-r)
3. Program description
(4) where b ~ - l ) is the element of the ith row and kth column of the matrix A r-1. The symbol s represents the number of equations in the system, the dimension of matrix A.
To minimize the sum of squares the program uses (with suitable modification according to our problem) the Gauss-Newton method given by Jennrich and Sampson in program BMDX85 [5]. Figure 1 gives the
For testing the accuracy of the Taylor-series method the error o f y may be approximated by:
where K is the number of steps in the integration of the system [7]. The advantage of the expression (5) is that only few calculations are required, but according to our experience it is a severe condition, therefore if it fails we use the method of halving the stepsize. The number of multiplications necessary to solve the system and to compute the ay/aakt partials can be
1
~'~[
1 STEP
[ I
PLOT
] I
Fig. 1. Structure of the program.
137
B. Kany~r, Z Erfdi, Parame~ of linear sys~ms without num~ical diffemnt~t~n PROBLEM IDENTIF.
ROSE-BENg./TIST
NuM§ER OF THE PARAMETERS
6
NUMBER OF THE OBSERVATIONs
19
NUMOER OF THE iNDEP,VAR|ABLES
1
NUMBER OF THE ~EPENCE NT VARIAB,
2
HuMB|R OF THE dEIGHT
|
VARIAB,
VE|GHT CODE
"1 *1
NUMBER OF THE O I F F . E ~ * °S
4
NUMBER OF TkE ?RIRER PAR,oS
4
1NP~T OATS FORMAT
|
I
(ISFS,I)
NUMBER OF MAXIMUM FUNCTION EVAL,
21
ToLEB,OF FUNCTIDM CHANG|NG
D,|OIEBe
T~LER,CR|TER, AT ~ATR*INvERS,
Q,BIII|I
N~MBER OF THE TAYLOR POWER
4
TA START,VALUE OF IkOEP,VAR
Q,P
ERRCR LIMIT AT soL,OF DoE,S,
Bl*illl
STARTING VALUES oF THE STATE VAR.S INOEcES OF NON-ZERO NATR,ELES
I
lO|,l I 2
2 1
I*l ) 2
O,O
0.0
4 3
F~.2a. Problemdescription on output.
structure of the program with the following subroutines: MAIN The purpose of MAIN is to read input data and print the description of the problem (fig. 2a). Input data Number of the primary and secondary parameters, number of the differential equations, number of data points, number of curves to be fitted simultaneously, number of weight factors, weight codes in case of weighting by measured data, tolerance at nonlinear regression, initial values of the system, error limit in solution of the system, structure of the matrix A given by the indices of the nonzero matrix elements, data format, data, limits and guesses of the parameter values,
MATCO The MATCO subroutine computes the powers of matrix and the aAm/Oakl partials given in eq. (4).
SETDIF It solves the differential equation system and calculates the ay/aakt partial derivatives according to the formulas (2) and (3). Testing of the accuracy is made by relation (5) and by halving of stepsize. OBSEQ The observation equations are defined by supplied FORTRAN statements. The routine OBSEQ describes the connection of the model and response function, e.g., measuring the sum of two compartments.
B. Kany6r, .1..Er6di, Parameters o/linear systems without numerical differentiation
138
MINIMA
lEO
MAXIM~
l,idO|£
IT£RITIo~
ERROR M£AN S@uARE
0 O
7.6138E
O.g O0
l.Og|~E
0.0
O,O
L.OOegE OJ
l.Ji|OE
O]
03
3.d|OdE
03
5.~llE
ld
6.|RIOE 5.8839E
It El
1.0~0E-03 7,4635E-04
2.|~00E-|$ 2.3283E°10
2.0704E012 4.5612E°12
1.1462E002 2.)621E°E2
t,3543E-12 3.7336E-12
203283E-10 2,3283E-10
4.9904E | l 4,2928E O~ 4.5194[ Id
1.3275E-Et 6.6374EoR2 2.9~43E-02
6.0572E-02 6.9665E-02 l.geEtE-Ot
1.8668E002 l.i2d4E-|2 5,2714E-13
L,ES96E-02 1.7dlgE-|2 4,|2| E-E2
3.9SOLE 3,44|1E 2.O~64E
7.8SO6t-i2 ?.9567Eol2
2.1709E Oi
9.8740E OR
2,ill6E 2,OIllE
tO eO
9.9913E 9.9~05E
II O~
00 OI 0~
9.99|6E 9.99|6E 9.9906E 9.9916E 9.9906E 9.9906E
|d Ol J~ go El OI
l.S4S9E 6.3447E
7
O
2.3467E i |
2,4295E-E2
l.lO?9E-IL
7.|OIEE-O3
8 9
I |
3,4838E'13 9.8768E'14
2.5017E-12 2.5096E-O2
I.IIlIE-II l.ig16E-01
7.9699E-03 8*0262E-R3
II 1! 12 13 14
I ] 2
9.6507E'04 9.8506E*E4 9.8$1LE'O4
4
9.8529E°|4
2 4
9,8506E'04 9.8539E°04
2.5~91E-02 2.5096E-E2 2.5096E-|2 2.5g96E-R2 2.5~96E°02 2.5~96E-12
l,iil6E-II I.iII6E-iI I.OgI6E*I& I.I|L6E°O! t.ROI6EoOt 1.0|16£-0!
1S
02
5.illiE-g3 ~.1133E0|3
1 I
4 5 6
5.~OIOE
l.leOiE-02 1.|~85Eo02
| 1 1
2 )
OQ
0.|
PARAMETERS
7.5269E 03 I . t ~ 7 9 E 13 4 . 0 4 L 4 [ 02 3 . $ 6 6 9 E 02
1
O.O
12 It
7.8184E-R2
il 00 iO
6.R264E-03
7.6217E-R2
8.02~4E-03
?.0217E-R3
8.|2~4E-|3
7.8217E*02 ?.8217E-E2 7.8217E°|2 7.6217E-E2
2.1NllE 2.0Site 2.o~ltE 2.OOlIE 2.OItlE 2.lillE
T,4124Eo|4
1.6867E-03
6.0264E0|3 8.|264E°03 8.|264E-|3
oo IO OI
1.5259E005 6.2842E O~ 9.eII8E 9.7087E I.|151E
id O| O~
ASYMPTOTIC S~A~OARO OEVIATIONS OF THE PARAMETERS 4.9351EoE5
9.8151E-05
i,~483E-OS
4.41eSE-SJ
ASYMPTOTIC CORRELATIo N MATRIX oF THE PARAHETERS I
t 2
2
3
4
5
6
t.OO~il
0,79672
0*66579
00.77737
.0039473
°0.50638
3
0,79672 1.66879
l,liiJi e*417|1
0,41711 l,gliiO
od.6997g -~.667~4
-0,48t19 0E,06797
e0.72202 0~,43154
4 S
°0.77737 -i,39473
°|,067|4 °0.16797
0,10827
0.78647 *i,07?73
6
-i.SI638
*0.69970 -i*4Oll9 "1*72212
-i043154
L.Iilig 1018827 #078647
|,lifll -I,i7773
I,lilll
ESTIMATED MATRIX OF THE OIFF,EQeSYSTEM
o 8 . 1 1 J t 6 l • ! 2 S 1| d,|0~16 -|,0531~
~,i
•e
CASE
I. 0 I,O
l" ! O,O
e . e 78~ 2
o. 0
0. IIdSg3 °11,07822 e. |
RESIDUALS
Fig. 2b. The error me a n square and t he pa ra ma t e r values after
I,E
each iteration, t he standard deviations, correlations and the estimated parameter values after the last iteration.
iNPUT VAR|AGLES
53LUT|oN OF THE [IFF.EQ.SYSTEH
1 2 ] 4 5
i.|46 4.229 g.$62 8.3g9 0,49L
-i,176 -E.294 -|,528 5.223 0,043
l,ill 2.J16 ).14| 4,006 5.|00
916.111 623.481 7S|,glg 605.001 6200000
2752.w4| 3416.1Ja 4Beo,eoI 4520.cd8 4067.~d|
1 2 3 4 ~
l,l~l 2.Jd9 3.|8| 4,081 S,088
6 ? 6 9 18 tl 12 13 14
0.112 -|.|$2 -0.002 -8,260 -0.363 -d.672 -8.541 -8.411 0.273
-0,906 *i,227 "4.145 2.261 - e.688 2,?73 1,160 -2,|62 1,320
8.|06 10.000 15.000 28,J0| 380065 40.000 50.J8| 6~,|06 78,g00
493.480 42?.g80 318.000 257,006
6 7 e 9 le 11 12 12 14
8.880 10.|80 l$.Odo 20.odo
178,|0J 163.000 152.800 1430880
6016,~80 6StA.~O| 7276.~86 7614.e | 8 7644,~08 73660088 6986.~80 6584.~ 0| 6|96.~88
15 16 17 18 19
8.132 8.376 8.095 0.333 0.123
°1.778 -4.J27 4.270 L.633 -~,038
6JeLJO 900088 180*|0| 118.860 12~0696
134°100 126.088 118.000 111.888 194,060
5810.~86 5461.~8| 3135.~00 4821.~08 4316,~8
15 16 17 18 19
202.008
F ~ . 2e. Residuals and input data.
9~,588 U2.277 74.944 68.469 62.751 49,209
9.~76 17.576 24,736 3Q,969 36.436 48.860
1.93T i.z4e 8.294 0.488 0.712 1,493
J.~11 i.,j8 8,42d 0,~53 f,181 0.358
48.080 5J,JiJ 68.080 78,080
42.785 31.000 25.727 20.236 17.087 16.354 1~.241 14.273
54.611 63.129 66.587 66.676 63.852 60,319 560734 53.260
2.049 3.372 4.450 $~789 6.263 6.320 6.136 5.656
0.635 1.699 3.237 T.~96 12,~57 17.~07 21,~86 26.502
~.OdO 90.000 18~.88J 110.086 12~.868
13.387 12.562 L1,790 11,867 16.367
58.029 460962 44.062 41.377 36,036
5.542 50223 4.913 4.616 4,335
31.~41 35.251 39.215 42.940 46.489
3~.88!
F ~ . 2d. S ol ut i on of t he differential e q u a t i o n system a ~ e r the Mst iteration.
139
B. Kanydr, J. Erbdi, Parameters of linear systems without numerical differentiation PLOTS OF TME tEIIOUAL$ I*0109 2*Oil| 3*6J19 4°6166 5*1|9d 0o0990 19*|910 IS,OO|9 2 29*0g90
IO|$EmV,oPtIOICT.wg.) a
| a
)9*6099
114 YM! INO|P|UDtNT VaRIAOk| i | j | t ~ 1 I t I I 2
I
l
4|*1||9
1
5906990
2
I
60.6010
|
2
t
70,9600
I
6900999
~
90.0018
2
I
2
t
1J|,||90
I
119.6990
laO*60dO
|
2
t
2 I -5.d352
l *4.0094
1 -2.96)6
I
I
"1.9570
6 PUg~kF~|~ 4EO|ZTAS AtJ~VA I v l Z S I I N T E S N | I
*0.9320
2
t I
I
0.0937
1.1195
I 2.1453
I
t
).1711
4,1969
I S.22~7
0007d4 : t
Fig. 29. Plots of the residuaL. MINIZ and STEP These subroutines minimize the sum of squares, estimates the parameters, the standard deviations and the correlation matrix by use of the GaussNewton method. They are only slightly modified from the BMDX85 subroutines. MINIZ prints the values of the parameters and the solution of the differ°ntial equation system.
According to the output code the predicted values, the observed data and the residuals may be printed and plotted. The subroutine PLOT from IBM SSP library plots the defined data.
4. Sample run The compartment model used for this run is shown in fig. 3. The model is used to analyse the rose-bengal kinetics in the investigation of liver function [9]. That kind of model is used in our experimental investigations too. The tracer is injected into compartment 1. The rate equations of the model in matrix form are the following:
dy _i Or21
"--(O~12+ (x32)
0
cza2 0
--~43
~ -
•y [0
~x4a
B. Kanydr, J. Er6di, Parameters of linear systems without numerical differentiation
140 ( ~ - - ~ ( ~
"3___~2~
~,___.~3~
~,=I0 ~,,~~23 ~
/2~
f2
Fig. 3. Compartmental model used for testing the program. The three angles represent observations that are directly proportional to the contents of the compartments, the proportional coefficients are o t 1, o21, o22 and a 23. In case of our rose-bengal investigations two data sets fl(t i) and f2(ti) are to be fitted simultaneously. where y is the amount of tracer in the different compartments and akt are the transport coefficients, the primary parameters. The three angles represent observations and the secondary parameter oij is the fraction of compartment ] measured by detector i. To test the program simulated data sets (fl(ti) and f2(ti) ) were generated. The exact values of the data were rounded to three decimal digits. The values of the primary parameters were: al2 = 0.025
Or21 = 0.100
(~32 = 0.008
0t43 = 0.080
and the secondary parameters, the proportional coefficients o21 = 2 and 022 = 023 = 10. The proportional coefficient o1 z was considered as a constant value and not as a parameter, ozl = 10. Figures 2 a - e show the result of the run. The problem and model description is given in fig. 2a. The weight code (q) determines the weighting factors (wi) in the sum of squares, namely wi = f(ti) q. The next figure (2b) shows the error mean square and the parameter values after each iteration, the standard deviations, correlations and the estimated parameters in matrix form after the last iteration. The residuals and the input data, the variables time, fl(ti) and f2 (ti) are printed in fig. 2c. The program gives the solution of the differential equations (fig. 2d) and the plot of the residuals (fig. 2e), after the last iteration, According to the last iteration the estimated values of the parameters approach the original values quite well.
5. Hardware and software specifications
was tested on an IBM 360/75 system. In the case of a three-compartment model with four parameters [4] the total t processor h e gtime parameters. u wase 5 - 1s5 ssdepending e s on initial of the The new version of the program runs on an RIAD20 with DOS. Input unit is the card reader and the output is given at the printer. The required memory is 100 kbyte in case of 12 differential equations 15 parameters and 7 value of the power of Taylor-serie. DOS requires 20 kbyte. The sample run with 4 equations, 6 parameters and two data sets uses 10 min total time (read input, processing and print output).
6. Mode of availability The FORTRAN IV program list may be obtained from the authors at the following address: Semmelweis University of Medicine, Computer Group, 1089 Budapest, Kulich Gy. 5, Hungary.
References [1 ] Y. Bard, Nonlinear Parameter Estimation (Academic Press, New York, London, 1974). [2] M. Berman and M.F. Weiss, SAAM Manual (National Institutes of Health, Bethesda, 1971). [3] J. Buell, R. Kalaba, A. Yakush and E. Ruspini, Comput.
Prog. Biomed. 2 (1971) 8. [4] S. Gy~Srgyiand B. Kany~r, Acta Biochem. Biophys. (Acad. Sci. Hung) 7 (1972) 359. [5 ] R.I. Jennricb and P. Sampson, in: Biomedical Computer Programs (X-Series suppl.) ed. W.J. Dixon (Univ. Calif.
Press, 1972). [6] R.I. Jennrich and P.B. Bright, Technometrics 17 (1975) 245. [7] B. Kany~irand J. T6th, Alkalmazott Matematikai Lapok (Letters of Applied Mathematics) in Hungarian (1978) in press. [8] w. Murray, Numerical Methods for Unconstrained Optimization (Academic Press, New York, London, 1972). 19] A.D. Waxman, P.A. Leins and J.K. Siemsen, Comput. Biomed. Res. 5 (1973) 1.