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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \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]\displaystyle{ B }[/math], которая позволяла бы получить искомый вектор [math]\displaystyle{ X }[/math] с требуемыми характеристиками в результате линейного преобразования [math]\displaystyle{ X = BY + M }[/math], где [math]\displaystyle{ Y }[/math] — n-мерный случайный вектор с независимыми нормально распределенными компонентами со стандартными параметрами.

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

[math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ Y }[/math] независимы и имеют стандартные параметры, справедливо выражение:

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

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

[math]\displaystyle{ \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]\displaystyle{ \Sigma }[/math], а в правых — элементы искомой матрицы [math]\displaystyle{ B }[/math]. Последовательно решая эту систему, получаем формулы для расчета элементов [math]\displaystyle{ \sigma_{ij} }[/math]:

[math]\displaystyle{ 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]\displaystyle{ B }[/math] имеет вид:

[math]\displaystyle{ 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]

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

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

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

[math]\displaystyle{ 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]\displaystyle{ u_j^i - }[/math] реализации случайной величины [math]\displaystyle{ u }[/math], а [math]\displaystyle{ k - }[/math] параметр, обеспечивающий качество приближения.

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