' ******************************************************************
' Multinormal distribution
' ******************************************************************

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

DECLARE FUNCTION RanGaussStd AS DOUBLE
' Random number from the standard normal distribution

DECLARE FUNCTION RanGauss(Mu AS DOUBLE, Sigma AS DOUBLE) AS DOUBLE
' Random number from normal dist (mean Mu, SD Sigma)

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

SUB RanMult(M() AS DOUBLE, L() AS DOUBLE, X() AS DOUBLE)
' ------------------------------------------------------------------
' Generates a random vector X from a multinormal distribution.
' M is the mean vector, L is the Cholesky factor (lower triangular)
' of the variance-covariance matrix.
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub  ' Array bounds
  DIM AS INTEGER I, J    ' Loop variables

  Lb = LBOUND(M)
  Ub = UBOUND(M)

  DIM  AS DOUBLE U(Lb TO Ub)

  ' Form a vector U of independent standard normal variates
  FOR I = Lb TO Ub
    U(I) = RanGaussStd
  NEXT I

  ' Form X = M + L * U, which follows the multinormal distribution
  FOR I = Lb TO Ub
    X(I) = M(I)
    FOR J = Lb TO I
      X(I) = X(I) + L(I, J) * U(J)
    NEXT J
  NEXT I
END SUB

SUB RanMultIndep(M() AS DOUBLE, S() AS DOUBLE, X() AS DOUBLE)
' ------------------------------------------------------------------
' Generates a random vector X from a multinormal distribution with
' uncorrelated variables. M is the mean vector, S is the vector
' of standard deviations.
' ------------------------------------------------------------------

  DIM AS INTEGER I

  FOR I = LBOUND(M) TO UBOUND(M)
    X(I) = RanGauss(M(I), S(I))
  NEXT I
END SUB
