testmcmc.pas

Delphi 的数学控件

{ **********************************************************************
  *                         Program TESTMCMC.PAS                       *
  *                             Version 1.3d                           *
  *                     (c) J. Debord, February 2003                   *
  **********************************************************************
  This program simulates a multinormal distribution by Markov Chain
  Monte Carlo (MCMC) using the Hastings-Metropolis algorithm.

  Although MCMC is best used when there is no direct way to simulate
  the distribution, it is used here for demonstration purposes since
  its results can be compared to those of the direct method (program
  RANMUL.PAS).

  The pdf P(X) of the multinormal distribution is such that:

                        P(X) = C * Exp(- F(X) / T)

  where F(X) = (X - M)' * V^(-1) * (X - M)

        C = 1/sqrt(|V| * (2*Pi)^N)

        T = 2

        M is the mean vector and V the variance-covariance matrix of
        the distribution. N is the dimension of the distribution.

        The constant C is not used in the simulation.
        
  The mean vector and variance-covariance matrix are stored in a data
  file with the following structure:

    Line 1            : Title of study
    Line 2            : Number of variables (N), e.g. 2 for binormal
    Line 3 to (N + 2) : Means and standard deviations
    Next lines        : Correlation coefficients, in
                        lower triangular matrix form

  The file RANMUL.DAT is an example data file.

  The results are stored in the output file TESTMCMC.TXT
  ********************************************************************** }

program testmcmc;

uses
  fmath, matrices, mcmc;

const
  TEMP = 2.0;   { Temperature }

var
  Title   : String;   { Title of study }
  N       : Integer;  { Number of variables }
  M       : TVector;  { Mean vector of original distribution }
  V       : TMatrix;  { Variance-covariance matrix of original distribution }
  V_inv   : TMatrix;  { Inverse variance-covariance matrix }
  Xmat    : TMatrix;  { Matrix of simulated vectors }
  Msim    : TVector;  { Mean of simulated distribution }
  Vsim    : TMatrix;  { Variance-covariance matrix of simulated distrib. }
  X_min   : TVector;  { Coordinates of the minimum of F(X)
                        = mode of simulated distribution }
  F_min   : Float;    { Value of F(X) at minimum }
  I       : Integer;  { Loop variable }
  ErrCode : Integer;  { Error code }


  function ReadParam(FileName : String; var Title : String;
                     var N : Integer; var M : TVector;
                     var V, V_inv : TMatrix) : Integer;
  var
    F    : Text;     { Data file }
    I, J : Integer;  { Loop variables }
    S    : TVector;  { Standard deviations }
    R    : Float;    { Correlation coefficient }
  begin
    Assign(F, FileName);
    Reset(F);

    Readln(F, Title);
    Readln(F, N);

    DimVector(M, N);
    DimVector(S, N);
    DimMatrix(V, N, N);
    DimMatrix(V_inv, N, N);

    { Read means and standard deviations. Compute variances }
    for I := 1 to N do
      begin
        Read(F, M[I], S[I]);
        V[I][I] := Sqr(S[I]);
      end;

    { Read correlation coefficients and compute covariances }
    for I := 2 to N do
      for J := 1 to Pred(I) do
        begin
          Read(F, R);
          V[I][J] := R * S[I] * S[J];
          V[J][I] := V[I][J];
        end;

    { Compute the inverse of the variance-covariance matrix }
    ReadParam := InvMat(V, 1, N, V_inv);

    Close(F);
  end;

  function ObjFunc(X : TVector) : Float;
  { Computes the function F(X) }
  var
    Sum1, Sum2 : Float;
    I, J       : Integer;
    D          : TVector;
  begin
    DimVector(D, N);

    for I := 1 to N do
      D[I] := X[I] - M[I];

    Sum1 := 0.0;
    for I := 1 to N do
      Sum1 := Sum1 + V_inv[I][I] * Sqr(D[I]);

    Sum2 := 0.0;
    for I := 2 to N do
      for J := 1 to Pred(I) do
        Sum2 := Sum2 + V_inv[I][J] * D[I] * D[J];

    ObjFunc := Sum1 + 2.0 * Sum2;
  end;

  procedure WriteResults(Title : String; M : TVector;
                         V : TMatrix; N : Integer);
  var
    I, J : Integer;
    S    : TVector;
    R    : Float;
  begin
    WriteLn;
    WriteLn(Title);
    WriteLn;

    WriteLn('      Mean      S.D.');
    WriteLn('--------------------');

    DimVector(S, N);
    for I := 1 to N do
      begin
        S[I] := Sqrt(V[I][I]);
        Writeln(M[I]:10:4, S[I]:10:4);
      end;

    WriteLn;
    WriteLn('Correlation matrix:');
    WriteLn;

    for I := 2 to N do
      begin
        for J := 1 to Pred(I) do
          begin
            R := V[I][J] / (S[I] * S[J]);
            Write(R:10:4);
          end;
        WriteLn;
      end;
  end;

  procedure WriteOutputFile(Title : String; Xmat : TMatrix; N : Integer);
  var
    F    : Text;
    I, J : Integer;
  begin
    Assign(F, 'testmcmc.txt');
    Rewrite(F);

    WriteLn(F, Title);
    Write(F, ' Iter');
    for I := 1 to N do
      Write(F, '        X', I);
    WriteLn(F);

    for I := 1 to MH_SavedSim do
      begin
        Write(F, I:5);
        for J := 1 to N do
          Write(F, Xmat[J][I]:10:4);
        WriteLn(F);
      end;

    Close(F);
  end;

begin
  if ReadParam('ranmul.dat', Title, N, M, V, V_inv) = MAT_SINGUL then
    begin
      WriteLn('Variance-covariance matrix is singular!');
      Exit;
    end;

  DimVector(Msim, N);
  DimVector(X_min, N);
  DimMatrix(Vsim, N, N);
  DimMatrix(Xmat, N, MH_SavedSim);

  { Initialize the mean vector and the variance-covariance matrix.
    For the sake of demonstration we start at a distance from the
    true mean and with enhanced standard deviations. }
  for I := 1 to N do
    begin
      Msim[I] := 3.0 * M[I];
      Vsim[I,I] := 10.0 * V[I,I];
    end;

  { Perform Metropolis-Hastings simulations }
  Write('Running. Please wait...');
  ErrCode := Hastings(ObjFunc, TEMP, Msim, Vsim,
                      1, N, Xmat, X_min, F_min);
  if ErrCode = 0 then
    begin
      WriteResults('Original distribution', M, V, N);
      WriteResults('Simulated distribution', Msim, Vsim, N);
      WriteOutputFile(Title, Xmat, N);
    end
  else
    WriteLn('Variance-covariance matrix is not positive definite!');
end.