' ******************************************************************
' Monte-Carlo simulation of the statistical distribution
' of the regression parameters for 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)"

CONST FirstParam = 1   ' Index of first parameter
CONST LastParam  = 2   ' Index of last parameter

CONST Temp     = 2     ' Temperature
CONST NCycles  = 10    ' Number of cycles
CONST MaxSim   = 1000  ' Max nb of simulations
CONST SavedSim = 1000  ' Nb of saved simulations

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

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

' Parameters and bounds
DIM SHARED AS DOUBLE B(FirstParam TO LastParam)
DIM SHARED AS DOUBLE B_min(FirstParam TO LastParam)
DIM SHARED AS DOUBLE B_max(FirstParam TO LastParam)

' Optimized parameters
DIM AS DOUBLE B_opt(FirstParam TO LastParam)

' Variance-covariance matrix of parameters
DIM AS DOUBLE V(FirstParam TO LastParam, FirstParam TO LastParam)

' Simulated parameter values
DIM AS DOUBLE B_sim(1 TO SavedSim, FirstParam TO LastParam)

' Value of objective function at minimum
DIM AS DOUBLE F_opt

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

SUB InitParam(B() AS DOUBLE, B_min() AS DOUBLE, B_max() AS DOUBLE, V() AS DOUBLE)
' ------------------------------------------------------------------
' Initialize parameters, bounds and variance-covariance matrix
' ------------------------------------------------------------------

  DIM AS INTEGER I, J, Lb, Ub
  DIM AS DOUBLE  R

  Lb = LBOUND(B)
  Ub = UBOUND(B)

  B_min(1) = 100 : B_max(1) = 1000
  B_min(2) = 0.1 : B_max(2) = 1

  FOR I = Lb TO Ub
    R = B_max(I) - B_min(I)
    B(I) = B_min(I) + 0.5 * R
    FOR J = Lb TO Ub
      IF I = J THEN
        ' The parameter range is assumed to cover 6 SD's
        V(I,J) = R * R / 36
      ELSE
        V(I,J) = 0
      END IF
    NEXT J
  NEXT I
END SUB

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

  CONST MaxFunc = 1E+30

  DIM AS INTEGER I
  DIM AS DOUBLE  S

  ' Check if the parameters are inside the bounds
  FOR I = LBOUND(B) TO UBOUND(B)
    IF B(I) < B_min(I) OR B(I) > B_max(I) THEN
      ObjFunc = MaxFunc
      EXIT FUNCTION
    END IF
  NEXT I

  S = 0

  FOR I = 1 TO N

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

    R(I) = Y(I) - Ycalc(I)
    S = S + R(I) * R(I)

    IF S >= MaxFunc THEN
      ObjFunc = MaxFunc
      EXIT FUNCTION
    END IF

  NEXT I

  ObjFunc = S
END FUNCTION

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

SUB WriteResults(X() AS DOUBLE, Y() AS DOUBLE, Ycalc() AS DOUBLE, _
                 R() AS DOUBLE, B_opt() AS DOUBLE, V() AS DOUBLE, _
                 F_opt AS DOUBLE)
' ------------------------------------------------------------------
' Writes results to screen
' ------------------------------------------------------------------

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

  DIM AS INTEGER I, M
  DIM AS DOUBLE  SB, Sr

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

  PRINT "Parameter    Est.value         Std.dev."

  PRINT Line1

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

  PRINT Line1

  M = UBOUND(B) - LBOUND(B) + 1  ' Number of parameters
  Sr = SQR(F_opt / (N - M))      ' Residual standard deviation

  PRINT USING "Number of observations : n =####"; N
  PRINT USING "Residual error         : s =####.########"; Sr

  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(B_sim() AS DOUBLE)
' ------------------------------------------------------------------
' Plots the distribution of the regression parameters
' ------------------------------------------------------------------

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

  DIM AS DOUBLE X(1 TO SavedSim), Y(1 TO SavedSim)

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

  DIM AS INTEGER I, Lb, Lb1

  Lb = LBOUND(B_sim, 2)
  Lb1 = Lb + 1

  ' Extract coordinates
  FOR I = 1 TO SavedSim
    X(I) = B_sim(I, Lb)
    Y(I) = B_sim(I, Lb1)
  NEXT I

  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

  WriteTitles "Monte Carlo Simulation", "B" + STR$(Lb), "B" + STR$(Lb1)

  SetClipping

  COLOR 12

  PlotCurve X(), Y(), 2, 3

  CloseGraph
END SUB

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

' Initialize parameters and variance-covariance matrix

InitParam B(), B_min(), B_max(), V()

' Select RNG

' SetRNG RNG_MWC
SetRNG RNG_MT
' SetRNG RNG_UVAG

' Perform simulation

InitMHParams NCycles, MaxSim, SavedSim

PRINT "Running. Please wait..."

Hastings @ObjFunc, Temp, B(), V(), B_sim(), B_opt(), F_opt

IF MathErr <> MatOk THEN
  PRINT "Unable to perform simulation!"
  SLEEP
END IF

' Write results
WriteResults X(), Y(), Ycalc(), R(), B_opt(), V(), F_opt

' Plot graphics
PlotGraph B_sim()

SLEEP
