БЧА НИВЦ МГУ. GSN4R_C (GSN5R_C). Генерация случайных чисел

Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/g_htm_c/gsn4r_c.htm
Дата изменения: Wed Feb 25 17:28:08 2015
Дата индексирования: Sun Apr 10 00:52:01 2016
Кодировка: Windows-1251

Текст подпрограммы и версий
gsn4r_c.zip , gsn5r_c.zip

Тексты тестовых примеров
tgsn4r_c.zip , tgsn5r_c.zip

Подпрограмма: gsn4r_c (версия gsn5r_c)

Назначение

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

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

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

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

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

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

    int gsn4r_c (integer *iseed, integer *n, integer *k,
            real *sigma, real *rvec, real *wkvec, integer *ierr)

   Параметры

iseed -	целая переменная, значение которой перед обращением к подпрограмме может быть любым целым числом в пределах [1, 2147483646]; по окончании работы ей присваивается новое значение, котоpое может быть использовано при последующем вхождении в подпрограмму;
n -	заданное количество k - мерных псевдослучайных векто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 только левый нижний треугольник ковариационной матрицы). По окончании работы подпрограммы матрица ковариаций заменяется на матрицу L^T такую, что S = L*L^T;
rvec -	вещественный двумерный массив размерности n на k, содержащий полученные k - мерные псевдослучайные векторы;
wkvec -	вещественный вектоp длины k используемый подпрограммой для хранения промежуточных результатов;
ierr -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы, при этом:

ierr=65 -

когда ковариационная матрица SIGMA не является положительно определенной.

   Версии

gsn5r_c -

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

   Вызываемые подпрограммы

afh1r_c -	треугольное разложение положительно определенной симметричной матрицы методом квадратного корня (методом Холецкого) с компактной формой представления симметричной и треугольной матриц.
gsn1r_c -	генерация массива псевдослучайных чисел, ноpмально распределенных с нулевым средним значением и единичной дисперсией.
utgs10_c -	подпрограмма выдачи диагностических сообщений при работе подпрограмм gsn4r_c и gsn5r_c.

   Замечания по использованию

B подпрограмме gsn5r_c в качестве параметра sigma надо подавать не ковариационную матрицу, а матрицу L^T такую, что LL^T является ковариационной матрицей. Эта матрица получается в результате работы подпрограммы gsn4r_c. Поэтому, если нужно получить несколько массивов псевдослучайных чисел с заданной ковариационной матрицей, то первый раз должно быть обращение к gsn4r_c, а все последующие - к gsn5r_c.

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

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

int main(void)
{
    /* Initialized data */
    static float sigma[6] = { 1.f,0.f,1.f,0.f,0.f,1.f };

    /* Local variables */
    static float rvec[6] /* was [2][3] */;
    static int ierr;
    extern int gsn4r_c(int *, int *, int *, float *, float *, float *, int *);
    static int k, n, iseed, i;
    static float wkvec[3];

    iseed = 831670774;
    n = 2;
    k = 3;
    gsn4r_c(&iseed, &n, &k, sigma, rvec, wkvec, &ierr);

    for (i = 1; i <= 6; ++i) {
        printf("\n %16.7e \n",rvec[i-1]);
    }
    return 0;
} /* main */


Результаты:

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