' ******************************************************************
' Simulation by Markov Chain Monte Carlo (MCMC) with the
' Metropolis-Hastings algorithm.
'
' This algorithm simulates the probability density function (pdf)
' of a vector X. The pdf P(X) is written as:
'
'                      P(X) = C * Exp(- F(X) / T)
'
' Simulating P by the Metropolis-Hastings algorithm is equivalent to
' minimizing F by simulated annealing at the constant temperature T.
' The constant C is not used in the simulation.
'
' The series of random vectors generated during the annealing step
' constitutes a Markov chain which tends towards the pdf to be
' simulated.
'
' It is possible to run several cycles of the algorithm.
' The variance-covariance matrix of the simulated distribution is
' re-evaluated at the end of each cycle and used for the next cycle.
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE MatOk 0
' No error

#DEFINE MatNotPD -4
' Matrix is not positive definite

' ------------------------------------------------------------------
' Machine-dependent constant
' ------------------------------------------------------------------

#DEFINE MaxNum 1.797693134862315D+308
' Max. floating point number: 2^1024

#DEFINE MaxLog 709.782712893384#
' Max. argument for Exp

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

DECLARE SUB Cholesky(A() AS DOUBLE, L() AS DOUBLE)
' Cholesky factorization of a positive definite symmetric matrix

DECLARE FUNCTION RanGen3 AS DOUBLE
' Returns a random number in (0, 1)

DECLARE SUB RanMult(M() AS DOUBLE, L() AS DOUBLE, X() AS DOUBLE)
' Random vector X from multinormal dist. (mean M, Cholesky factor L)

DECLARE SUB RanMultIndep(M() AS DOUBLE, S() AS DOUBLE, X() AS DOUBLE)
' Random vector X from multinormal dist. (uncorrelated)

' ******************************************************************

DIM SHARED AS INTEGER MH_NCycles  = 10    ' Number of cycles
DIM SHARED AS INTEGER MH_MaxSim   = 1000  ' Max nb of simulations
DIM SHARED AS INTEGER MH_SavedSim = 1000  ' Nb of saved simulations

SUB InitMHParams(NCycles  AS INTEGER, _
                 MaxSim   AS INTEGER, _
                 SavedSim AS INTEGER)
' ------------------------------------------------------------------
' Initializes Metropolis-Hastings parameters
' ------------------------------------------------------------------

  IF NCycles  > 0 THEN MH_NCycles = NCycles
  IF MaxSim   > 0 THEN MH_MaxSim = MaxSim
  IF SavedSim > 0 AND SavedSim <= MaxSim THEN MH_SavedSim = SavedSim
END SUB

SUB CalcSD(V() AS DOUBLE, S() AS DOUBLE)
' ------------------------------------------------------------------
' Computes the standard deviations for independent random vectors
' from their variance-covariance matrix.
' ------------------------------------------------------------------

  DIM AS INTEGER I, Lb, Ub

  Lb = LBOUND(V, 1)
  Ub = UBOUND(V, 1)

  I = Lb
  ErrCode = MatOk

  DO
    IF V(I,I) > 0 THEN
      S(I) = SQR(V(I,I))
    ELSE
      ErrCode = MatNotPD
    END IF
    I = I + 1
  LOOP UNTIL (ErrCode <> MatOk) or (I > Ub)
END SUB

FUNCTION Accept(DeltaF AS DOUBLE, T AS DOUBLE) AS INTEGER
' ------------------------------------------------------------------
' Checks if a variation DeltaF of the function at temperature T is
' acceptable.
' ------------------------------------------------------------------

  DIM AS DOUBLE X

  IF DeltaF < 0.0 THEN
    Accept = -1
    EXIT FUNCTION
  END IF

  X = DeltaF / T

  IF X > MaxLog THEN  ' Exp(-X) ~ 0
    Accept = 0
    EXIT FUNCTION
  END IF

  Accept = (EXP(- X) > RND)
END FUNCTION

SUB HastingsCycle(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                        T      AS DOUBLE, _
                        X()    AS DOUBLE, _
                        V()    AS DOUBLE, _
                        Indep  AS INTEGER, _
                        Xmat() AS DOUBLE, _
                        Xmin() AS DOUBLE, _
                  BYREF Fmin   AS DOUBLE)
' ------------------------------------------------------------------
' Performs one cycle of the Metropolis-Hastings algorithm
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub         ' Array bounds
  DIM AS INTEGER I, J, K        ' Loop variables
  DIM AS INTEGER Iter           ' Iteration count
  DIM AS INTEGER FirstSavedSim  ' Index of 1st simul. to be saved
  DIM AS DOUBLE  F, F1          ' Function values
  DIM AS DOUBLE  DeltaF         ' Variation of function
  DIM AS DOUBLE  Sum            ' Statistical sum

  Lb = LBOUND(X)
  Ub = UBOUND(X)

  DIM AS DOUBLE S(Lb TO Ub), X1(Lb TO Ub), L(Lb TO Ub, Lb TO Ub)

  ' Compute SD's or Cholesky factor from var-cov matrix
  IF Indep THEN
    CalcSD V(), S()
  ELSE
    Cholesky V(), L()
  END IF

  IF ErrCode = MatNotPD THEN EXIT SUB

  ' Compute initial function value
  F = Func(X())

  ' Perform MH_MaxSim simulations at constant temperature
  FirstSavedSim = MH_MaxSim - MH_SavedSim + 1
  Iter = 1
  K = 1

  DO
    ' Generate new vector
    IF Indep THEN
      RanMultIndep X(), S(), X1()
    ELSE
      RanMult X(), L(), X1()
    END IF

    ' Compute new function value
    F1 = Func(X1())
    DeltaF = F1 - F

    ' Check for acceptance
    IF Accept(DeltaF, T) THEN
      PRINT ".";
      FOR I = Lb TO Ub
        X(I) = X1(I)
      NEXT I

      IF Iter >= FirstSavedSim THEN
        ' Save simulated vector into line K of matrix Xmat
        FOR I = Lb TO Ub
          Xmat(K, I) = X1(I)
        NEXT I
        K = K + 1
      END IF

      IF F1 < Fmin THEN
        ' Update minimum
        FOR I = Lb TO Ub
          Xmin(I) = X1(I)
        NEXT I
        Fmin = F1
      END IF

      F = F1
      Iter = Iter + 1
    END IF
  LOOP UNTIL Iter > MH_MaxSim

  ' Update mean vector
  FOR I = Lb TO Ub
    Sum = 0
    FOR K = 1 TO MH_SavedSim
      Sum = Sum + Xmat(K,I)
    NEXT K
    X(I) = Sum / MH_SavedSim
  NEXT I

  ' Update variance-covariance matrix
  FOR I = Lb TO Ub
    FOR J = I TO Ub
      Sum = 0
      FOR K = 1 TO MH_SavedSim
        Sum = Sum + (Xmat(K,I) - X(I)) * (Xmat(K,J) - X(J))
      NEXT K
      V(I,J) = Sum / MH_SavedSim
    NEXT J
  NEXT I

  FOR I = Lb + 1 TO Ub
    FOR J = Lb TO I - 1
      V(I,J) = V(J,I)
    NEXT J
  NEXT I
END SUB

SUB Hastings(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                   T      AS DOUBLE, _
                   X()    AS DOUBLE, _
                   V()    AS DOUBLE, _
                   Xmat() AS DOUBLE, _
                   Xmin() AS DOUBLE, _
             BYREF Fmin   AS DOUBLE)
' ------------------------------------------------------------------
' Simulation of a probability density function by the
' Metropolis-Hastings algorithm
' ------------------------------------------------------------------
' Input parameters :  Func = Function such that the pdf is
'                                P(X) = C * Exp(- Func(X) / T)
'                     T    = Temperature
'                     X    = Initial mean vector
'                     V    = Initial variance-covariance matrix
' ------------------------------------------------------------------
' Output parameters : Xmat = Matrix of simulated vectors, stored
'                            by lines, i.e.
'                            Xmat(1..MH_SavedSim, Lb..Ub)
'                     X    = Mean of distribution
'                     V    = Variance-covariance matrix of distrib.
'                     Xmin = Coordinates of minimum of Func(X)
'                               (mode of the distribution)
'                     Fmin = Value of Func(X) at minimum
' ------------------------------------------------------------------
' The function MathErr() returns one of the following codes:
'   MatOk    = no error
'   MatNotPD = var-cov matrix is not positive definite
' ------------------------------------------------------------------

  DIM AS INTEGER K, Indep

  ' Initialize the random number generator
  RANDOMIZE TIMER

  K = 1
  Indep = -1
  Fmin = MaxNum

  DO
    HastingsCycle Func, T, X(), V(), Indep, Xmat(), Xmin(), Fmin
    Indep = 0
    K = K + 1
  LOOP UNTIL ErrCode <> MatOk or K > MH_NCycles
END SUB

