Участник:Lonalone/Генерация гауссовского вектора методом линейных преобразований

Автор описания: Меньших И. М.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Из многомерных распределений особый интерес представляет нормальное. Этому закону подчиняются все результаты воздействия большого числа случайных факторов, среди которых нет превалирующих.

В статье приведен алгоритм генерации n-мерного гауссовского случайного вектора с помощью метода линейных преобразований. Известно, что случае нормально распределенного случайного вектора, ковариационная матрица вместе с математическим ожиданием этого вектора полностью определяют его распределение. Поэтому для полного статистического соответствия моделируемого и теоретического распределения гауссовского вектора достаточно обеспечить требуемые значения указанных параметров.

Идея алгоритма заключается в линейном преобразовании n-мерного случайного вектора [math]Y[/math], компоненты которого независимы и одинаково распределены по нормальному закону со стандартными параметрами, в случайный вектор [math]X[/math] с требуемыми ковариационной матрицей и вектором математических ожиданий.

1.2 Математическое описание алгоритма

1.2.1 Метод линейных преобразований

Даны ковариационная матрица [math]\Sigma[/math] и вектор математических ожиданий [math]M[/math]:

[math] \Sigma = \|\sigma_{ij}\| = \| \mathbb{E}[(X_{i} - m_{x_{i}})(X_{j} - m_{x_{j}})]\|, \\ M = (m_{x_{1}}, m_{x_{2}}, ..., m_{x_{n}})^T. [/math]

Требуется найти такую матрицу [math]B[/math], которая позволяла бы получить искомый вектор [math]X[/math] с требуемыми характеристиками в результате линейного преобразования [math]X = BY + M[/math], где [math]Y[/math] — n-мерный случайный вектор с независимыми нормально распределенными компонентами со стандартными параметрами.

Будем искать матрицу [math]B[/math] в виде нижней треугольной матрицы. Перейдем от матричной записи к системе алгебраических уравнений:

[math] \begin{pmatrix} X_{1} \\ X_{2} \\ \vdots \\ X_{n} \end{pmatrix} = \begin{pmatrix} b_{11} & 0 & \cdots & 0 \\b_{21} & b_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ b_{n1} & b_{n2} & \cdots & b_{nn} \end{pmatrix} \times \begin{pmatrix} Y_{1} \\ Y_{2} \\ \vdots \\ Y_{n} \end{pmatrix} + \begin{pmatrix} m_{x_{1}} \\ m_{x_{2}} \\ \vdots \\ m_{x_{n}} \end{pmatrix} \Rightarrow [/math]

[math] \begin{cases}X_{1} - m_{x_{1}} = b_{11}Y_{1} \\X_{2} - m_{x_{2}} = b_{21}Y_{1} + b_{22}Y_{2} \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\X_{n} - m_{x_{n}} = b_{n1}Y_{1} + b_{n2}Y_{2} + \cdots + b_{nn}Y_{n} \end{cases} [/math]

Поскольку компоненты вектора [math]Y[/math] независимы и имеют стандартные параметры, справедливо выражение:

[math] \mathbb{E}[Y_{i}Y_{j}] = \left\{\begin{matrix} 1, &i = j, \\ 0, &i \not= j. \end{matrix}\right. [/math]

Почленно перемножив сами на себя и между собой соответственно левые и правые части уравнений системы и взяв от результатов перемножения математическое ожидание, получаем систему уравнений вида:

[math] \begin{cases} \mathbb{E}[(X_{1} - m_{x_{1}})(X_{1} - m_{x_{1}})] = \mathbb{E}[b_{11}Y_{1}b_{11}Y_{1}], \\ \mathbb{E}[(X_{1} - m_{x_{1}})(X_{2} - m_{x_{2}})] = \mathbb{E}[(b_{21}Y_{1} + b_{22}Y_{2})b_{11}Y_{1}], \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \end{cases} [/math]

Таким образом, в левых частях полученной системы уравнений располагаются элементы заданной ковариационной матрицы [math]\Sigma[/math], а в правых — элементы искомой матрицы [math]B[/math]. Последовательно решая эту систему, получаем формулы для расчета элементов [math]\sigma_{ij}[/math]:

[math] b_{11}=\sqrt{\sigma_{11}}; b_{21}=\frac{\sigma_{12}}{\sqrt{\sigma_{11}}}; b_{22} = \sqrt{\sigma_{22} - \frac{\sigma_{12}}{\sigma_{11}}}, \cdots [/math]

Рекуррентная формула для расчета любого элемента матрицы преобразования [math]B[/math] имеет вид:

[math] b_{ij} = \frac {\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}} {\sqrt{\sigma_{ij} - \sum_{k=1}^{j-1} b_{jk}^2}}, \quad 1 \leqslant j \leqslant i \leqslant n [/math]

(суммы с верхним нулевым пределом считаются равными нулю)

Таким образом, матрица [math]B[/math] получается с помощью разложения Холецкого матрицы [math]\Sigma[/math].

1.2.2 Генерация случайного вектора Y

Пусть имеется генератор псевдослучайных чисел (ГПСЧ, PRNG), с помощью которого можно получить реализацию случайной величины [math]u \sim U(0,1) [/math]. Описанный выше случайный вектор [math]Y=(y_1, \cdots, y_n)[/math] с независимыми компонентами составим из [math]n[/math] реализаций случайной величины [math]\eta \sim N(0,1)[/math]. Каждую такую реализацию [math]y_i[/math], в свою очередь, получим с помощью приближения по ЦПТ [math]k[/math] случайными величинами, распределенными равномерно на отрезке [0,1]:

[math] y_i = \frac {1} {\sqrt{\frac{k}{12}} } (\sum_{j=1}^k{u_j^i} - \frac{k}{2}), \forall i = 1, \cdots, n \\ [/math]

где [math]u_j^i - [/math] реализации случайной величины [math]u[/math], а [math]k - [/math] параметр, обеспечивающий качество приближения.

Как правило, [math]k[/math] берут равным 12 и считают, что для подавляющего числа практических задач обеспечивается должная точность вычислений.

1.3 Вычислительное ядро алгоритма

Вычисление выражений:
[math]\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}, \quad 1 \leqslant j \leqslant i \leqslant n[/math]

1.4 Макроструктура алгоритма

• Заполнение матрицы [math]B[/math]
• Генерация вектора [math]Y[/math]
• Вычисление вектора [math]X = BY + M[/math]

1.5 Схема реализации последовательного алгоритма

1) Заполнение матрицы [math]B[/math].

(a) [math]b_{11}= \sqrt{\sigma_{11}}[/math]

(b) [math]b_{k1}= \frac{\sigma_{k1}}{l_{11}}[/math], при [math]k = \overline{2,n}[/math]

Следующие два пункта выполняются циклически, друг за другом для [math]i = \overline{2,n}[/math] (переход к новому [math]i[/math] после пункта 4).

(c) [math]b_{ii} = \sqrt{\sigma_{ii} - \sum_{p = 1}^{i - 1} b_{ip}^2}[/math]

(d) (кроме [math]i = n[/math]): [math]b_{ji} = \frac{\sigma_{ji} - \sum_{p = 1}^{i - 1} b_{ip} b_{jp}} {l_{ii}}[/math], при [math]j = \overline{i+1,n}[/math]

причем вычисление числителя дроби осуществляется в режиме накопления.

2) Генерация вектора [math]Y[/math].

Следующие два пункта выполняются циклически, друг за другом для [math]i = \overline{1,n}[/math].

(a) [math]u_j^i \leftarrow PRNG,[/math] при [math]j = \overline{1,k}[/math]

(b) [math]y_i = \frac {1} {\sqrt{\frac{k}{12}} } (\sum_{j=1}^k{u_j^i} - \frac{k}{2})[/math]

3) Вычисление вектора [math]X[/math].

[math]x_i = m_{x_i} + \sum_{k = 1}^{i} b_{ik} y_{k},[/math] при [math]i = \overline{1,n}[/math].

1.9 Входные и выходные данные алгоритма

• Входные данные: вещественная всюду плотная положительно определенная симметрическая [math](n \times n)[/math] матрица [math]\Sigma[/math] и вещественный [math](n)[/math] вектор [math]M[/math];
• Выходные данные: вещественный [math](n)[/math] вектор [math]X[/math].

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

Вариант последовательной реализации без режима накопления:

// compute matrix B (core)
for (int i = 0; i < n; i++){
    for (int j = 0; j <= i; j++){
        double sum_up = 0, sum_down = 0;
        for (int k = 0; k < j; k++){
            sum_up += matrix_B[i][k] * matrix_B[j][k];
            sum_down += matrix_B[j][k] * matrix_B[j][k];
        }
        matrix_B[i][j] = (matrix_cov[i][j] - sum_up) / sqrt(matrix_cov[j][j] - sum_down);
    }
}
// generate vector Y
vec_Y = generate_Y();
// linear transformation
for (int i = 0; i < n; i++){  
    vec_X[i] = vec_M[i];
    for (int k = 0; k <= i; k++){
        vec_X[i] += matrix_B[i][k] * vec_Y[k];
    }
}

Вариант последовательной реализации с режимом накопления:

// compute matrix B (core)
matrix_B = low_triangle(matrix_cov);
for (int i = 0; i < n; i++){
    matrix_B[i][i] = sqrt(matrix_B[i][i]);
    for (int j = i+1; j < n; j++){
        matrix_B[j][i] /= matrix_B[i][i];
    }
    // accumulate
    for (int k = i+1; k < n; k++){
        for (int j = k; j < n; j++){
            matrix_B[j][k] -= matrix_B[j][i] * matrix_B[k][i];
        }
    }
}
// generate vector Y
vec_Y = generate_Y();
// linear transformation
for (int i = 0; i < n; i++){  
    vec_X[i] = vec_M[i];
    for (int k = 0; k <= i; k++){
        vec_X[i] += matrix_B[i][k] * vec_Y[k];
    }
}

2.7 Существующие реализации алгоритма

Несмотря на то, что точечный метод Холецкого реализован в множестве пакетов (напр., LINPACK, LAPACK, SCALAPACK), моделирование гауссовского случайного вектора не имеет такой популярности. MS Excel и распространённые статистические пакеты (напр., SPSS, Statistica) позволяют моделировать только одномерные статистические распределения.

Таким образом, имеется возможность составить многомерное распределение из нескольких одномерных, но только при условии, что компоненты независимы. Если же нужно исследовать данные с зависящими друг от друга переменными, приходится работать над написанием программы.