*-----------------------------------------------------------------------
* Fit Ratkowsky's pasture growth data
*-----------------------------------------------------------------------
      PROGRAM PASTUR
       IMPLICIT NONE
       INTEGER N
       PARAMETER (N = 9)
       INTEGER J, IERR
       DOUBLE PRECISION A, B, C, B1, B2, B3,
     $            B1EST, B2EST, B3EST, DET, RELERR
       DOUBLE PRECISION X(N), Y(N), ODY(N), W0(N), W1(N), Z(N), 
     $                   ZV(N), ZG(N), ZK(N)
       DOUBLE PRECISION DN, RESULT, SX, SX2, SXY, SY, T
       DOUBLE PRECISION ZERO, ONE
       PARAMETER (ZERO = 0.D0, ONE = 1.D0)

*-----------------------------------------------------------------------
* NIST/ITL StRD
* Dataset Name:  Rat42             (Rat42.dat)
*
* Procedure:     Nonlinear Least Squares Regression
*
* Description:   This model and data are an example of fitting
*                sigmoidal growth curves taken from Ratkowsky (1983).
*                The response variable is pasture yield, and the
*                predictor variable is growing time.
*
*
* Reference:     Ratkowsky, D.A. (1983).
*                Nonlinear Regression Modeling.
*                New York, NY:  Marcel Dekker, pp. 61 and 88.
*
* Data:          1 Response  (y = pasture yield)
*                1 Predictor (x = growing time)
*                9 Observations
*                Higher Level of Difficulty
*                Observed Data
*
* Model:         Exponential Class
*                3 Parameters (b1 to b3)
*
*                y = b1 / (1+exp[b2-b3*x])  +  e
*
*
*           Starting Values                  Certified Values
*
*         Start 1     Start 2           Parameter     Standard Deviation
*   b1 =   100         75            7.2462237576E+01  1.7340283401E+00
*   b2 =     1          2.5          2.6180768402E+00  8.8295217536E-02
*   b3 =     0.1        0.07         6.7359200066E-02  3.4465663377E-03
*
* Residual Sum of Squares:                    8.0565229338E+00
* Residual Standard Deviation:                1.1587725499E+00
* Degrees of Freedom:                                6
* Number of Observations:                            9
*-----------------------------------------------------------------------
       DATA Y(1), X(1) /   8.930E0,        9.000E0 /
       DATA Y(2), X(2) /  10.800E0,       14.000E0 /
       DATA Y(3), X(3) /  18.590E0,       21.000E0 /
       DATA Y(4), X(4) /  22.330E0,       28.000E0 /
       DATA Y(5), X(5) /  39.350E0,       42.000E0 /
       DATA Y(6), X(6) /  56.110E0,       57.000E0 /
       DATA Y(7), X(7) /  61.730E0,       63.000E0 /
       DATA Y(8), X(8) /  64.620E0,       70.000E0 /
       DATA Y(9), X(9) /  67.080E0,       79.000E0 /

       DATA B1 / 7.2462237576E+01 /
       DATA B2 / 2.6180768402E+00 /
       DATA B3 / 6.7359200066E-02 /

*-----------------------------------------------------------------------
* Fit the transformed model
*
*      1/y = 1/b1 + exp(b2)/b1 exp(-b3 x)
*
* Thus:
*         z:  1/y
*         a:  1/b1
*         b:  exp(b2)/b1
*         c:  -b3
*
*     z = a + b exp(c x)
*
*-----------------------------------------------------------------------
  20   FORMAT ()
  21   FORMAT ('...', A, '...')
       PRINT 21, 'Logistic model'
       DO J = 1, N
         ODY(J) = ONE / Y(J)
       END DO
       CALL RATEST (X, ODY, N, C, W1, IERR)
       PRINT *, 'ratest returns #', IERR, ' estimated c: ', C

       SX = ZERO                  ! least squares line
       SY = ZERO
       SX2 = ZERO
       SXY = ZERO
       DO J = 1, N
         T = EXP(C * X(J))
         SX = SX + T
         SY = SY + ODY(J)
         SX2 = SX2 + T**2
         SXY = SXY + T * ODY(J)
       END DO
       DN = DBLE(N)

       DET = DN * SX2 - SX**2
       A = (SY * SX2 - SX * SXY) / DET
       B = (N * SXY - SX * SY) / DET

       B1EST = 1 / A
       RELERR = (B1EST - B1) / B1
       PRINT *, 'estimated b1: ', B1EST, ' certified b1: ', B1,
     $    ' relative error: ', RELERR

       B2EST = LOG(B1EST * B)
       RELERR = (B2EST - B2) / B2
       PRINT *, 'estimated b2: ', B2EST, ' certified b2: ', B2,
     $    ' relative error: ', RELERR

       B3EST = -C
       RELERR = (B3EST - B3) / B3
       PRINT *, 'estimated b3: ', B3EST, ' certified b3: ', B3,
     $   ' relative error: ', RELERR

*-----------------------------------------------------------------------
*             Error of the estimate
*-----------------------------------------------------------------------
       DO J = 1, N                        ! integral method
         ZV(J) = B1EST / (ONE + EXP(B2EST - B3EST * X(J)))
       END DO
       CALL CODD (N, Y, ZV, RESULT)
       PRINT *, 'Integral method coefficient of determination: ', RESULT

       DO J = 1, N
         Z(J) = B1 / (ONE + EXP(B2 - B3 * X(J)))
       END DO
       CALL CODD (N, Y, Z, RESULT)
       PRINT *, 'coefficient of determination, certified parameters: ',
     $               RESULT

*-----------------------------------------------------------------------
*             Try Gompertz model.  Transform:
*
*      y = a exp(-b exp(-c x))
*      ln(y) = ln(a) - b exp(-c x)
*
* Thus:
*        z:  ln(y)
*        a':  ln(a)
*        b':  -b
*        c':  -c
*
*      z = a + b exp(c x)
*
*-----------------------------------------------------------------------
       PRINT 20
       PRINT 21, 'Gompertz model'
       DO J = 1, N
         ODY(J) = LOG(Y(J))
       END DO
       CALL RATEST (X, ODY, N, C, W1, IERR)
       PRINT *, 'ratest returns #', IERR, ' parameter c: ', C

       SX = ZERO                  ! least squares line
       SY = ZERO
       SX2 = ZERO
       SXY = ZERO
       DO J = 1, N
         T = EXP(C * X(J))
         SX = SX + T
         SY = SY + ODY(J)
         SX2 = SX2 + T**2
         SXY = SXY + T * ODY(J)
       END DO
       DN = DBLE(N)

       DET = DN * SX2 - SX**2
       A = (SY * SX2 - SX * SXY) / DET
       B = (N * SXY - SX * SY) / DET

       A = EXP(A)
       B = -B
       C = -C
       PRINT *, 'a: ', A,  ' b: ', B

       DO J = 1, N
         ZG(J) = A * EXP(-B * EXP(-C * X(J)))
       END DO
       CALL CODD (N, Y, ZG, RESULT)
       PRINT *, 'coefficient of determination, Gompertz model: ', RESULT

*-----------------------------------------------------------------------
*             Try Korf model.
*
*      y = a exp(-b x^(-c))
*
* Subroutine KOREST internally applies the transformation
*      ln y = (ln a) - b x^(-c)
*      z:  ln(y)
*      a':  ln(a)
*      c':  c
* given the original (X,Y) values, and returns growth rate C.
*
*-----------------------------------------------------------------------
       PRINT 20
       PRINT 21, 'Korf model'
       CALL KOREST (X, Y, N, C, W0, W1, IERR)
       PRINT *, 'korest returns #', IERR, ' parameter c: ', C

       SX = ZERO                  ! least squares line
       SY = ZERO
       SX2 = ZERO
       SXY = ZERO
       DO J = 1, N
         T = X(J)**(-C)
         SX = SX + T
         SY = SY + ODY(J)
         SX2 = SX2 + T**2
         SXY = SXY + T * ODY(J)
       END DO
       DN = DBLE(N)

       DET = DN * SX2 - SX**2
       A = (SY * SX2 - SX * SXY) / DET
       B = (N * SXY - SX * SY) / DET

       A = EXP(A)
       B = -B
       PRINT *, 'a: ', A,  ' b: ', B

       DO J = 1, N
         ZK(J) = A * EXP( -B * X(J)**(-C) )
       END DO
       CALL CODD (N, Y, ZK, RESULT)
       PRINT *, 'coefficient of determination, Korf model: ', RESULT

       PRINT *, 'x, y, z_cert, Verhulst, Gompertz, Korf'
       DO J = 1, N
         PRINT *, X(J), Y(J), Z(J), ZV(J), ZG(J), ZK(J)
       END DO

       PRINT *, 'program complete'
      END  ! of pastur