C-------------------------------- DEMO ----------------------------------------- C C Note: the original code source has been lost, this one was typed again by C JRL on 03-DEC-2006 from an hard copy of the article: C C @article{LobryJR1991, C Author = {Lobry, J.R. and Rosso, L. and Flandrois, J.-P.}, C Title = {A {FORTRAN} subroutine for the determination of parameter C confidence limits in non-linear models}, C Journal = {Binary}, C Volume = {3}, C Pages = {86-93}, C Year = {1991} } C C This program compiles for me without warning with g77 : C GNU Fortran (GCC) 3.4.2 C C VARIABLES C---------- C C SSR NAME OF AN USER SUPPLIED SUBROUTINE WHICH RETURNS THE SUM OF C SQUARED RESIDUALS BETWEEN OBSERVED AND THEORETICAL VALUES C C PN THE NUMBER OF PARAMETERS IN THE MODEL C C POPT VECTOR WHICH CONTAINS THE THE PARAMETER VALUES SUCH THAT C THE SSR IS MINIMUM C C PINF VECTOR WHICH CONTAINS THE PARAMETER LOWER BOUNDS TO C GENERATE STARTING POINTS C C PSUP VECTOR WHICH CONTAINS THE PARAMETER UPPER BOUNDS TO C GENERATE STARTING POINTS C C THRES THRESHOL VALUE FOR THE SSR WHICH MAY BE OBTAINED FROM: C C THRES = SSR(POPT)*(1+PN*F(alpha;PN,N-PN)/(N-PN)), C C WHERE N IS THE NUMBER OF POINTS IN THE DATA SET. C ACCORDING TO BEALE (1960): C C @article{BealeEML1960, C Author = {Beale, E.M.L.}, C Title = {Confidence regions in non-linear estimation}, C Journal = {Journal of the Royal Statistical Society}, C Year = {1960}, C Volume = {22B}, C Pages = {41-88} } C C MAXPT MAXIMUM NUMBER OF POINTS TO BE GENERATED C C SEED AN INITIAL CONDITION FOR THE PSEUDO-RANDOM NUMBERS C GENERATOR C C FILN THE NAME OF THE OUTPUT FILE TO WRITE COORDINATES OF POINTS C C EPS DETERMINES THE ACCURACY OF RESULTS DURING CONVERGENCE TO THE C THRESHOLD, THE STOP CRITERION IS THAT (THRES-SSR)/THRES < EPS C C MAXSSR A STARTING POINT IS REJECTED WHEN THE RATIO OF THE SSR AT THIS C POINT TO THE THRESHOLD VALUE IS GREATER THAN MAXSSR C (SSR/THRES > MAXSSR). THIS IS TO AVOID STARTING FROM A POINT TOO C FAR FROM THE OPTIMAL POINT: IN THIS CASE CONVERGENCE IS TOO C LONG AND SUBJECT TO FAILURE DUE TO NUMERICAL IMPRECISION C C MAXIT MAXIMUM NUMBER OF ITERATIONS ALLOWED FOR CONVERGENCE TO C THRESHOLD C C------------------------------------------------------------------------------- PROGRAM DEMO EXTERNAL SSR REAL*8 SSR,POPT(8),PINF(8),PSUP(8),THRES,SEED,EPS,MAXSSR INTEGER PN,MAXPT,MAXIT CHARACTER*30 FILN C C INITIALISATION OF VARIABLES C PN=4 C POPT(1): MuOpt PINF(1)=1.3 POPT(1)=1.3915 PSUP(1)=1.5 C POPT(2): Tmin PINF(2)=285.5 POPT(2)=289.398 PSUP(2)=292.0 C POPT(3): Topt PINF(3)=311.8 POPT(3)=313.25 PSUP(3)=314.5 C POPT(4): Tmax PINF(4)=319.9 POPT(4)=320.234 PSUP(4)=320.8 THRES=0.0395811*(1+(4.0*3.36)/(15.0-4.0)) MAXPT=10000 SEED=0.1071966 FILN="OUTDEMO" EPS=0.0001 MAXSSR=1.0E+05 MAXIT=500 C C CALLING SUBROUTINE CRUNCH: C CALL CRUNCH(SSR,PN,POPT,PINF,PSUP,THRES, & MAXPT,SEED,FILN,EPS,MAXSSR,MAXIT) END C----------------------------------- SSR -------------------------------------- C C PURPOSE C ------- C THIS FUNCTION RETURNS THE SUM OF SQUARED RESIDUALS BETWEEN OBSERVED C AND THEORETICAL VALUES. FOR CLARITY, EXPERIMENTAL DATA HAVE BEEN C INCLUDED DIRECTLY INSIDE THE PROGRAM. IT IS MORE USUAL, HOWEVER, TO C READ DATA FROM A FILE AT THE LEVEL OF THE MAIN PROGRAM AND THEN TO C PASS VALUES TO THE SUBROUTINE SSR VIA COMMON VARIABLES. C C Note added in Dec 2006: the included dataset of specific growth rate C versus temperature for Escherichia coli was deduced from figure 2 C in the paper: C C @article{RatkowskyDA1983, C Author = {Ratkowsky, D.A. and Lowry, R.K. and McMeekin, T.A. and C Stokes, A.N. and Chandler, R.E.}, C Title = {Model for bacterial culture growth rate throughout the entire C biokinetic range}, C Journal = {Journal of Bacteriology}, C Volume = {154}, C Pages = {1222-1226}, C Year = {1983} } C C COMPUTATION OF THE SSR IN THE CASE OF THE CARDINAL TEMPERATURE MODEL C T1 = Muopt, T2 = Tmin, T3 = Topt, T4 = Tmax C C VARIABLES C --------- C C T INPUT VECTOR WHICH CONTAINS PARAMETER VALUES C C P INPUT VARIABLE WHICH CONTAINS THE NUMBER OF PARAMETERS IN THE MODEL C C------------------------------------------------------------------------------- REAL*8 FUNCTION SSR(T,P) INTEGER P REAL*8 T(8) REAL*8 A,B,T1,T2,T3,T4 T1=T(1) T2=T(2) T3=T(3) T4=T(4) A=T2+T4 B=T2*T4 SSR= &(0.25-T1*(1-((294-T3)**2)/((294-T3)**2+294*(A-294)-B)))**2+ &(0.56-T1*(1-((296-T3)**2)/((296-T3)**2+296*(A-296)-B)))**2+ &(0.61-T1*(1-((298-T3)**2)/((298-T3)**2+298*(A-298)-B)))**2+ &(0.79-T1*(1-((300-T3)**2)/((300-T3)**2+300*(A-300)-B)))**2+ &(0.94-T1*(1-((302-T3)**2)/((302-T3)**2+302*(A-302)-B)))**2+ &(1.04-T1*(1-((304-T3)**2)/((304-T3)**2+304*(A-304)-B)))**2+ &(1.16-T1*(1-((306-T3)**2)/((306-T3)**2+306*(A-306)-B)))**2+ &(1.23-T1*(1-((308-T3)**2)/((308-T3)**2+308*(A-308)-B)))**2+ &(1.36-T1*(1-((310-T3)**2)/((310-T3)**2+310*(A-310)-B)))**2+ &(1.32-T1*(1-((312-T3)**2)/((312-T3)**2+312*(A-312)-B)))**2+ &(1.36-T1*(1-((314-T3)**2)/((314-T3)**2+314*(A-314)-B)))**2+ &(1.34-T1*(1-((316-T3)**2)/((316-T3)**2+316*(A-316)-B)))**2+ &(0.96-T1*(1-((318-T3)**2)/((318-T3)**2+318*(A-318)-B)))**2+ &(0.83-T1*(1-((319-T3)**2)/((319-T3)**2+319*(A-319)-B)))**2+ &(0.16-T1*(1-((320-T3)**2)/((320-T3)**2+320*(A-320)-B)))**2 END C--------------------------- CRUNCH ------------------------------------------- C C PURPOSE C ------ C THE SUBROUTINE CRUNCH WRITES COORDINATES OF POINTS SUCH THAT C THE SSR FUNCTION EQUALS THE THRESHOLD VALUE IN THE OUTPUT FILE C C------------------------------------------------------------------------------ SUBROUTINE CRUNCH(SSR,PN,POPT,PINF,PSUP,THRES, & MAXPT,SEED,FILN,EPS,MAXSSR,MAXIT) EXTERNAL SSR REAL*8 SSR,POPT(8),PINF(8),PSUP(8),THRES,SEED,EPS,MAXSSR INTEGER PN,MAXPT,MAXIT CHARACTER*30 FILN C C DECLARATION OF LOCAL VARIABLES C REAL*8 COOR(8),TETA(8),DPINF(8),DPSUP(8),COORSV(8),TETASV(8) REAL*8 LOWER,UPPER,SSRLOW,SSRUP,NEW,SSRNEW INTEGER ITERNBR LOGICAL OK C C CHECKING POPT VALUES C IF(SSR(POPT,PN).GT.THRES)THEN WRITE(*,*)"ERROR FROM ROUTINE CRUNCH: SSR IS GREATER" WRITE(*,*)"THAN THRESHOLD WITH OPTIMAL PARAMETER VALUES" WRITE(*,*)"THRESHOLD = ",THRES WRITE(*,*)"SSR = ",SSR(POPT,PN) STOP END IF C C OPENING OUTPUT FILE C OPEN(UNIT=34,FILE=FILN,STATUS="NEW") C C INTERMEDIARY COMPUTATIONS C DO 50,I=1,PN DPINF(I)=(POPT(I)-PINF(I))*SQRT(DBLE(PN)) DPSUP(I)=(PSUP(I)-POPT(I))*SQRT(DBLE(PN)) 50 CONTINUE C C MAIN LOOP OF THE PROGRAM. C REPETITION OF THE SAME WITH MAXPT DIFFERENT STARTING POINTS C I=1 1 I=I+1 C C CHOOSING A POINT AT RANDOM AT THE SURFACE OF A HYPERSHERE C OF RADIUS 1.0 WITH ITS CENTER AT THE ORIGIN: C CALL SPHERE(PN,COOR,SEED) C C COMPUTATION OF THE CORRESPONDING PARAMETER VALUES (SCALING) C DO 110,J=1,PN IF(COOR(J).GT.0.0)THEN TETA(J)=POPT(J)+COOR(J)*DPSUP(J) ELSE TETA(J)=POPT(J)+COOR(J)*DPSUP(J) END IF 110 CONTINUE SSRUP=SSR(TETA,PN) OK=.TRUE. C C TESTING IF THE POINT IS OUTSIDE THE CONFIDENCE REGION: C IF(SSRUP.LT.THRES)THEN WRITE(*,*)"WARNING 1 FROM SUBROUTINE CRUNCH" WRITE(*,*)"SSR LES THAN THRESHOLD WITH:" WRITE(*,*)"PARAM VALUE % DISTANCE FROM OPTIMAL PARAM VALUE" DO 120,K=1,PN WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN)) 120 CONTINUE WRITE(*,*)"-----------------------------------------" OK=.FALSE. I=I-1 END IF C C TESTING IF SSR/THRES IS NOT TOO BIG: C IF(SSRUP/THRES.GT.MAXSSR)THEN WRITE(*,*)"WARNING 2 FROM SUBROUTINE CRUNCH" WRITE(*,*)"SSR/THRES GREATER THAN:",MAXSSR WRITE(*,*)"PARAM VALUE % DISTANCE FROM OPTIMAL PARAM VALUE" DO 125,K=1,PN WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN)) 125 CONTINUE WRITE(*,*)"-----------------------------------------" OK=.FALSE. I=I-1 END IF C C SAVING COORDINATES C DO 115,K=1,PN TETASV(K)=TETA(K) COORSV(K)=COOR(K) 115 CONTINUE C C CONVERGENCE TO THRESHOLD, INITIALIZATION: C IF(OK)THEN LOWER=0.0 UPPER=1.0 DO 130,J=1,PN TETA(J)=POPT(J) 130 CONTINUE SSRLOW=SSR(TETA,PN) ITERNBR=0 C C CONVERGENCE TO THRESHOLD, MAIN LOOP: C 200 NEW=LOWER+(UPPER-LOWER)*SQRT((THRES-SSRLOW)/(SSRUP-SSRLOW)) ITERNBR=ITERNBR+1 DO 150,J=1,PN IF(COOR(J).GT.0.0)THEN TETA(J)=POPT(J)+NEW*COOR(J)*DPSUP(J) ELSE TETA(J)=POPT(J)+NEW*COOR(J)*DPINF(J) END IF 150 CONTINUE SSRNEW=SSR(TETA,PN) IF(SSRNEW.GT.THRES)THEN UPPER=NEW ELSE LOWER=NEW END IF IF(ITERNBR.LT.MAXIT)THEN IF(ABS(THRES-SSRNEW)/THRES.LT.EPS)THEN WRITE(34,1000)(TETA(K),K=1,PN) ELSE GO TO 200 END IF ELSE WRITE(*,*)"WARNING 3 FROM SUBROUTINE CRUNCH" WRITE(*,*)"NUMBER OF ITERATIONS GREATER THAN",MAXIT WRITE(*,*)"PARAM VALUE % DISTANCE FROM OPTIMAL PARAM VALUE" DO 135,K=1,PN WRITE(*,1000)TETA(K),100*COOR(K)*SQRT(DBLE(PN)) 135 CONTINUE WRITE(*,*)"-----------------------------------------" I=I-1 END IF END IF IF(I.LE.MAXPT)GO TO 1 CLOSE(34) 1000 FORMAT('',8(F12.6)) C C Note added in Dec. 2006: I'm not sure that this output format is C correct, I can't read well it on the hard copy of the listing. C END C------------------------------ SPHERE ---------------------------------------- C C PURPOSE C ------- C THE SUBROUTINE SPHERE RETURNS RANDOM COORDINATES FOR A POINT C BELONGING TO A HYPERSPHERE OF RADIUS 1.0 WITH ITS CENTER C AT THE ORIGIN (0,0,...,0) C C VARIABLES C --------- C N SPACE DIMENSION: NUMBER OF COORDINATES TO GENERATE. C N IS IN THE RANGE [2,8] C C COOR COORDINATES OF POINTS RETURNED BY SPHETEST (Note added in C Dec. 2006: this should be a typo for SPHERE) C C SEED SEED FOR THE PSEUDO-RANDOM GENERATOR, C ITS VALUE HAS TO BE INITIALIZED BY THE CALLING PROGRAM C WITHIN [0-1] AND NOT FURTHER MODIFIED C C COMMENTS C -------- C THE N COORDINATES OF A POINT ARE CHOSEN AT RANDOM WITHIN A C HYPERCUBE. THE DISTANCE FROM THIS POINT TO THE ORIGIN (0,0,...,0) IS C COMPUTED. IF THE DISTANCE IS LESS THAN 1.0, THEN THE POINT IS INSIDE THE C HYPERSHERE, AND THE COORDINATES ARE NORMALIZED SO THAT THE C DISTANCE TO THE ORIGIN EQUALS ONE. WHEN THE POINT IS OUTSIDE THE C HYPERSPHERE, ANOTHER POINT IS CHOSEN AT RANDOM. C THIS IS TO OBTAIN A HOMOGENEOUS DISTRIBUTION OF THE C POINTS ON THE SURFACE OF THE HYPERSPHERE. C C THE OBSERVED MEAN NUMBER OF TRIALS TO OBTAIN A POINT WITHIN A C HYPERSHERE IS GIVEN HERE AS FUNCTION OF THE SPACE DIMENSION N: C C N NUMBER OF TRIALS (MEAN OVER 1000) C 2 1.284 C 3 1.970 C 4 3.243 C 5 6.010 C 6 12.193 C 7 26.699 C 8 62.023 C C------------------------------------------------------------------------------- SUBROUTINE SPHERE(N,COOR,SEED) INTEGER N REAL*8 COOR(8),SEED,SUM 5 SUM=0.0 C C CHOICE OF N COORDINATES AT RANDOM C DO 10,I=1,N COOR(I)=2.0*RANDOM(SEED)-1.0 SUM=SUM+COOR(I)*COOR(I) 10 CONTINUE C C TESTING WHETHER THE POINT IS INSIDE THE HYPERSHERE: C IF(ABS(SQRT(SUM)).LE.1.0)THEN C C Note added in Dec. 2006: I don't understand why the absolute value is C taken in the line right above. Extreme defensive programing? A sum C of squares is always positive, this shouldn't be useful. Anyway, it C doesn't hurt, I keep it as in the original code. C DO 20,I=1,N COOR(I)=COOR(I)/SQRT(SUM) 20 CONTINUE ELSE GO TO 5 END IF END C------------------------------ RANDOM ---------------------------------------- C C PURPOSE C ------- C THE FUNCTION RANDOM RETURNS A PSEUDO-RANDOM NUMBER WITHIN [0-1] C C VARIABLES C --------- C C X SEED FOR THE PSEUDO-RANDOM GENERATOR, C ITS VALUE HAS TO BE INITIALIZED BY THE CALLING PROGRAM C WITHIN [0-1] AND NOT FURTHER MODIFIED C C NOTE ADDED IN DEC 2006 C ---------------------- C C This generator is known to behave satisfactorily for up to 100,000 C successive calls: C C @article{ChasseJL1974, C Author = {Chass{\'e}, J.-L. and Debouzie, D.}, C Title = {Utilisation des tests de {K}iveliovitch et {V}ialar dans C l'{\'e}tude de quelques g{\'e}n{\'e}rateurs de nombres pseudo- C al{\'e}atoires}, C Journal = {Revue de Statistique Appliqu{\'e}e}, C Volume = {22}, C Pages = {83-90}, C Year = {1974} } C C Enter ?".Random.seed" at your R prompt for references to better C pseudo-random generators. C C------------------------------------------------------------------------------ FUNCTION RANDOM(X) REAL*8 X,RANDOM RANDOM=MOD(262144*3125*X,262144.0)/262144.0 X=RANDOM END