' ******************************************************************
' This program performs a nonlinear least squares fit of the
' exponential model: Y = B(1) * Exp(- B(2) * X)
' ******************************************************************

#INCLUDE "fbmath.bi"
#INCLUDE "fbgfx.bi"

' ******************************************************************
' Special declarations (modify as necessary)
' ******************************************************************

' Number of observations
CONST N = 10

' Data
DIM SHARED AS DOUBLE X(1 TO N) = {  1,   2,   3,   4,   5,   6,  7,  8,  9, 10}
DIM SHARED AS DOUBLE Y(1 TO N) = {416, 319, 244, 188, 144, 113, 85, 66, 50, 41}

' Function name
CONST FuncName = "Y = B(1) * Exp(- B(2) * X)"

' Indexes of first and last fitted parameters
CONST FirstParam = 1, LastParam = 2

' Max number of iterations
CONST MaxIter = 1000

' Required precision
CONST Tol = 1E-6

' Significance level
CONST Alpha = 0.05

' ******************************************************************
' General declarations
' ******************************************************************

' Computed Y values and residuals
DIM SHARED AS DOUBLE Ycalc(1 TO N)
DIM SHARED AS DOUBLE R(1 TO N)

' Fitted parameters
DIM SHARED AS DOUBLE B(FirstParam TO LastParam)

' Gradient and Variance-covariance matrix
DIM SHARED AS DOUBLE G(FirstParam TO LastParam)
DIM SHARED AS DOUBLE V(FirstParam TO LastParam, FirstParam TO LastParam)

DIM AS DOUBLE Fmin    ' Value of objective function at minimum
DIM AS DOUBLE Det     ' Determinant of Hessian at minimum
DIM AS DOUBLE T, F    ' Student's t and Snedecor's F

DIM AS TRegTest Test  ' Statistical tests

' ******************************************************************
' Special subprograms (modify as necessary)
' ******************************************************************

SUB ApproxFit
' ------------------------------------------------------------------
' Approximate fit of the exponential model by linear regression:
' Ln(Y) = Ln(B(1)) - B(2) * X
' ------------------------------------------------------------------

  DIM AS INTEGER P                      ' Nb. of points for linear reg.
  DIM AS INTEGER K                      ' Loop variable

  P = 0                                 ' Count the number of points
  FOR K = 1 TO N                        ' which can be transformed
    IF Y(K) > 0 THEN P = P + 1
  NEXT K

  DIM AS DOUBLE X1(1 TO P), Y1(1 TO P)  ' Transformed coordinates
  DIM AS DOUBLE S1(1 TO P)              ' Standard deviations

  P = 0
  FOR K = 1 TO N
    IF Y(K) > 0 THEN
      P = P + 1
      X1(P) = X(K)
      Y1(P) = LOG(Y(K))
      S1(P) = 1 / Y(K)
    END IF
  NEXT K

  WLinFit X1(), Y1(), S1(), B(), V()

  IF MathErr = MatOk THEN
    B(1) = EXP(B(1))
    B(2) = - B(2)
  END IF
END SUB

FUNCTION ObjFunc(B() AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Objective function to be minimized
' ------------------------------------------------------------------

  DIM AS INTEGER K
  DIM AS DOUBLE  S

  S = 0

  FOR K = 1 TO N

    ' Regression function (change if necessary)
    Ycalc(K) = B(1) * EXP(- B(2) * X(K))

    R(K) = Y(K) - Ycalc(K)
    S = S + R(K) * R(K)
  NEXT K

  ObjFunc = S
END FUNCTION

SUB HessGrad(B() AS DOUBLE, G() AS DOUBLE, H() AS DOUBLE)
' ------------------------------------------------------------------
' Approximate Hessian and Gradient
' ------------------------------------------------------------------

  DIM AS INTEGER I, J, K

  DIM AS DOUBLE D(FirstParam TO LastParam)

  FOR I = FirstParam TO LastParam
    G(I) = 0
    FOR J = I TO LastParam
      H(I,J) = 0
    NEXT J
  NEXT I

  FOR K = 1 TO N

    ' Derivatives (change if necessary)
    D(1) = EXP(- B(2) * X(K))
    D(2) = - B(1) * X(K) * D(1)

    FOR I = FirstParam TO LastParam
      G(I) = G(I) - D(I) * R(K)
      FOR J = I TO LastParam
        H(I,J) = H(I,J) + D(I) * D(J)
      NEXT J
    NEXT I
  NEXT K

  FOR I = FirstParam + 1 TO LastParam
    FOR J = FirstParam TO I - 1
      H(I,J) = H(J,I)
    NEXT J
  NEXT I
END SUB

FUNCTION PltFunc(X AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Function to be plotted
' ------------------------------------------------------------------

  PltFunc = B(1) * EXP(- B(2) * X)
END FUNCTION

' ******************************************************************
' General subprograms
' ******************************************************************

SUB WriteResults(X() AS DOUBLE, Y() AS DOUBLE, Ycalc() AS DOUBLE, _
                 R() AS DOUBLE, B() AS DOUBLE, V() AS DOUBLE, _
                 Test AS TRegTest, T AS DOUBLE, F AS DOUBLE)
' ------------------------------------------------------------------
' Writes results to screen
' ------------------------------------------------------------------

  DIM AS STRING Line1 = STRING$(74, "-")
  DIM AS STRING Line2 = STRING$(74, "=")

  DIM AS INTEGER I
  DIM AS DOUBLE  SB, Sr

  PRINT Line2
  PRINT "Nonlinear regression: "; FuncName
  PRINT Line1

  PRINT "Parameter   Est.value         Std.dev.         ";
  PRINT USING "##"; 100 * (1 - Alpha);
  PRINT "% Confidence Interval"

  PRINT Line1

  FOR I = FirstParam TO LastParam
    SB = SQR(V(I,I))
    PRINT USING "B(#)"; I;
    PRINT USING "########.########"; B(I); SB; B(I) - T * SB;
    PRINT ";";
    PRINT USING "########.########"; B(I) + T * SB
  NEXT I

  PRINT Line1

  Sr = SQR(Test.Vr)

  PRINT USING "Number of observations           : n           = ###"; N
  PRINT USING "Residual error                   : s           = ###.########"; Sr
  PRINT USING "Coefficient of correlation       : r           = ###.########"; SQR(Test.R2)
  PRINT USING "Coefficient of determination     : r^2         = ###.########"; Test.R2
  PRINT USING "Adjusted coef. of determination  : r^2(a)      = ###.########"; Test.R2a
  PRINT USING "Variance ratio (explained/resid) : F(###, ###) = #######.####"; Test.Nu1; Test.Nu2; Test.F
  PRINT USING "Critical variance ratio          : F(p = #.##) = #######.####"; 1 - Alpha; F

  PRINT Line1

  PRINT "  i        Y obs.       Y calc.      Residual      Std.res."

  PRINT Line1

  FOR I = 1 TO N
    R(I) = Y(I) - Ycalc(I)
    PRINT USING "###"; I;
    PRINT USING "#########.####"; Y(I); Ycalc(I); R(I); R(I) / Sr
  NEXT I

  PRINT Line2
END SUB

SUB PlotGraph(X() AS DOUBLE, Y() AS DOUBLE, B() AS DOUBLE)
' ------------------------------------------------------------------
' Plots observed points and fitted curve
' ------------------------------------------------------------------

  DIM AS DOUBLE Xmin, Xmax, Xstep  ' Ox scale
  DIM AS DOUBLE Ymin, Ymax, Ystep  ' Oy scale

  IF NOT InitScreen THEN
    PRINT "Unable to set graphic mode!"
    EXIT SUB
  END IF

  InitGraph

  AutoScale X(), LinScale, Xmin, Xmax, Xstep
  AutoScale Y(), LinScale, Ymin, Ymax, Ystep

  SetOxScale LinScale, Xmin, Xmax, Xstep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 4
  PlotOyAxis Xmin, 4

  PlotGrid

  WriteTitles "Nonlinear regression: " + FuncName, "X", "Y"

  SetClipping

  COLOR 12

  PlotCurve X(), Y(), 1, 3

  PlotFunc @PltFunc

  CloseGraph
END SUB

' ******************************************************************
' Main program
' ******************************************************************

' Perform regression

ApproxFit

Marquardt @ObjFunc, @HessGrad, B(), MaxIter, Tol, Fmin, G(), V(), Det

IF MathErr <> MatOk THEN
  PRINT "Error in function minimization!"
  SLEEP
END IF

' Update variance-covariance matrix and compute statistical tests
RegTest Y(), Ycalc(), V(), Test

' Compute Student's t and Snedecor's F
T = InvStudent(Test.Nu2, 1 - 0.5 * Alpha)
F = InvSnedecor(Test.Nu1, Test.Nu2, 1 - Alpha)

' Write results
WriteResults X(), Y(), Ycalc(), R(), B(), V(), Test, T, F

' Plot graphics
PlotGraph X(), Y(), B()
