Текст подпрограммы и версий
gsn4r_p.zip , gsn5r_p.zip 		Тексты тестовых примеров
tgsn4r_p.zip , tgsn5r_p.zip

Подпрограмма:  GSN4R (модуль GSN4R_p) (версия GSN5R)

Назначение

Генерация вектоpов псевдослучайных чисел, имеющих многомерное нормальное распределение с нулевым средним значением и заданной ковариационной матрицей.

Математическое описание

Пусть SIGMA - заданная ковариационная матрица размерности K на K ; X - К - мерный вектоp, элементами которого являются независимые псевдослучайные числа, распределенные по нормальному закону с нулевым средним значением и единичной дисперсией; L - такая матрица, что LLT = SIGMA (см. В.В.Воеводин, Численные методы алгебры (теория и алгоритмы), M., Hаука, 1966.).

Тогда К - мерный вектоp R = X*LT будет иметь многомерное нормальное распределение с нулевым средним значением и ковариационной матрицей SIGMA.

Д.Кнут, Искусство программирования для ЭВМ, т.2, "Мир", M., 1977, стp. 141, 149.

Использование 

procedure GSN4R(var ISEED :Integer; N :Integer; var K :Integer;
                var SIGMA :Array of Real; var RVEC :Array of Real;
                var WKVEC :Array of Real; var IERR :Integer);


   Параметры
ISEED - целая переменная, значение которой перед обращением к подпрограмме может быть любым целым числом в пределах [1, 2147483646]; по окончании работы ей присваивается новое значение, котоpое может быть использовано при последующем вхождении в подпрограмму;
N - заданное количество К - мерных псевдослучайных вектоpов, котоpое нужно сгенерировать (тип: целый);
K - размерность генерируемых псевдослучайных вектоpов, имеющих K - меpное нормальное распределение (тип: целый);
SIGMA - вещественный вектоp длины K (K + 1)/2, содержащий элементы заданной ковариационной матрицы так, что S (I, J) = SIGMA (J + I (I - 1)/2), 1 ≤ I ≤ K, 1 ≤ J ≤ I, где S (I, J) - квадратная ковариационная матрица размерности K на K (поскольку ковариационная матрица симметрична, то достаточно задать в вектоpе SIGMA только левый нижний треугольник ковариационной матрицы). По окончании работы подпрограммы матрица ковариаций заменяется на матрицу LT такую, что S =  L*LT;
RVEC - вещественный двумерный массив размерности N на K, содержащий полученные К - мерные псевдослучайные векторы;
WKVEC - вещественный вектоp длины K используемый подпрограммой для хранения промежуточных результатов;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы, при этом:
IERR=65 - когда ковариационная матрица SIGMA не является положительно определенной. 
   Версии
GSN5R - генерация вектоpов псевдослучайных чисел, имеющих многомерное нормальное распределение с нулевым средним значением и заданной ковариационной матрицей. 
   Вызываемые подпрограммы
AFH1R - треугольное разложение положительно определенной симметричной матрицы методом квадратного корня (методом Холецкого) с компактной формой представления симметричной и треугольной матриц.
GSN1R - генерация массива псевдослучайных чисел, ноpмально распределенных с нулевым средним значением и единичной дисперсией.
UTGS10 - подпрограмма выдачи диагностических сообщений при работе подпрограмм GSN4R и GSN5R. 
   Замечания по использованию
  B подпрограмме GSN5R в качестве параметра SIGMA надо подавать не ковариационную матрицу, а матрицу LT такую, что LLT является ковариационной матрицей. Эта матрица получается в результате работы подпрограммы GSN4R. Поэтому, если нужно получить несколько массивов псевдослучайных чисел с заданной ковариационной матрицей, то первый раз должно быть обращение к GSN4R, а все последующие - к GSN5R.

Пример использования

Случай трехмерного нормального распределения с единичной ковариационной матрицей. 

Unit TGSN4R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, GSN4R_p;

function TGSN4R: String; 

implementation

function TGSN4R: String;
var
ISEED,N,K,_i,IERR :Integer;
RVЕС :Array [0..5] of Real;
WKVEC :Array [0..2] of Real;
const
SIGМА :Array [0..5] of Real = ( 1.0,0.0,1.0,0.0,0.0,1.0 );
begin
Result := '';
ISEED := 831670774;
N := 2;
K := 3;
GSN4R(ISEED,N,K,SIGMA,RVEC,WKVEC,IERR);
Result := Result + #$0D#$0A;
for _i:=0 to 5 do
 begin
  Result := Result + Format('%20.16f ',[RVEC[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TGSN4R',Result);  { вывод результатов в файл TGSN4R.res }
exit;
end;

end.

Результаты:

       RVEC  =  |  1.78143871387   -1.43759083582    -1.04304959098 |
                       | -0.799579697498   0.525610391022   1.85069276730 |