ranmull.pas

Delphi 的数学控件

{ **********************************************************************
  *                         Program RANMULL.PAS                        *
  *                            Version 1.4d                            *
  *                     (c) J. Debord, February 2003                   *
  **********************************************************************
  This program simulates a multi-lognormal distribution. The mean vector
  and variance-covariance matrix are stored in a data file with the
  following structure:

    Line 1            : Name of distribution.
    Line 2            : Size of distribution (N), e.g. 2 for bi-lognormal
    Line 3 to (N + 2) : Means and standard deviations.
    Next lines        : Correlation coefficients, in
                        lower triangular matrix form.

  The file RANMULL.DAT is an example data file.

  The results are stored in an output file (one random vector by line).
  ********************************************************************** }

program ranmull;

uses
  fmath, randnum, matrices, stat;

const
  NSIM = 400;  { Number of simulations }

var
  Name : String;   { Name of distribution }
  N    : Integer;  { Dimension of distribution }
  M    : TVector;  { Mean vector of original lognormal distribution }
  V    : TMatrix;  { Variance-covariance matrix of original lognormal dist. }
  M0   : TVector;  { Mean vector of auxiliary normal dist. }
  V0   : TMatrix;  { Variance-covariance matrix of auxiliary normal dist. }
  L    : TMatrix;  { Cholesky factor of V0 }
  Z    : TVector;  { Random vector from the auxiliary normal dist. }
  X    : TVector;  { Random vector from the original lognormal dist. }
  F    : Text;     { Output file }
  I, J : Integer;  { Loop variables }

  procedure ReadParam(FileName : String; var Name : String; var N : Integer;
                      var M : TVector; var V : TMatrix);
  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, Name);
    Readln(F, N);

    DimVector(M, N);
    DimVector(S, N);
    DimMatrix(V, 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;

    Close(F);
  end;

begin
  { Read parameters of log-normal distribution LN(M, V) }
  ReadParam('ranmull.dat', Name, N, M, V);

  DimVector(X, N);
  DimVector(Z, N);
  DimVector(M0, N);
  DimMatrix(V0, N, N);
  DimMatrix(L, N, N);

  { Define auxiliary normal distribution N(M0, V0) }
  for I := 1 to N do
    begin
      for J := 1 to N do
        V0[I,J] := Ln(V[I,J] / (M[I] * M[J]) + 1.0);
      M0[I] := Ln(M[I]) - 0.5 * V0[I,I];
    end;

  Assign(F, 'ranmull.out');
  Rewrite(F);

  { Perform Cholesky decomposition of variance-covariance matrix }
  if Cholesky(V0, 1, N, L) = MAT_NOT_PD then
    begin
      WriteLn('Variance-covariance matrix is not positive definite.');
      Exit;
    end;

  for I := 1 to NSIM do
    begin
      { Pick random vector from auxiliary normal distribution }
      RanMult(M0, L, 1, N, Z);

      { Convert to lognormal }
      for J := 1 to N do
        X[J] := Exp(Z[J]);

      { Output result to file }
      for J := 1 to N do
        Write(F, X[J]:12:6);
      Writeln(F);
    end;
  Close(F);
end.