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

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

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

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

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

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

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

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

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

$\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.$

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

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

$\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$

$\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}$

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

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

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

$\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}$

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

$b_{11}=\sqrt{\sigma_{11}}; b_{21}=\frac{\sigma_{12}}{\sqrt{\sigma_{11}}}; b_{22} = \sqrt{\sigma_{22} - \frac{\sigma_{12}}{\sigma_{11}}}, \cdots$

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

$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$

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

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

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

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

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

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

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

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

Вычисление следующих выражений (в режиме накопления или без него):
$\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}, \quad 1 \leqslant j \leqslant i \leqslant n$

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

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

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

1) Заполнение матрицы $B$.

(a) $b_{11}= \sqrt{\sigma_{11}}$

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

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

(c) $b_{ii} = \sqrt{\sigma_{ii} - \sum_{p = 1}^{i - 1} b_{ip}^2}$

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

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

2) Генерация вектора $Y$.

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

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

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

3) Вычисление вектора $X$.

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

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

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

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 и др.