Приложение С (справочное). Формирование выборок из распределений вероятностей

Приложение С
(справочное)

Формирование выборок из распределений вероятностей

С.1 Общие положения

С.1.1 В настоящем приложении приведены рекомендации по формированию выборки в соответствии с заданной функцией распределения вероятностей. Формирование такой выборки представляет собой ключевой момент при трансформировании распределений с использованием метода Монте-Карло. В качестве источников информации можно использовать [37] (сборник таблиц математических функций) и [38] (библиотека соответствующих программ).

С.1.2 Генератор псевдослучайных чисел для любого распределения, в том числе, для рассмотренных в 6.4 (см. таблицу 1), может теоретически быть получен на основе заданной функции распределения и генератора выборки для равномерного распределения (см. С.2). Генератор для равномерного распределения рассматривается в С.3.3. Для некоторых распределений, таких как нормальное распределение или t-распределение, более эффективным будет использование генераторов, специально разработанных именно для этих распределений (например, рассматриваемых в настоящем приложении). Общие рекомендации по формированию выборки в соответствии с заданным законом распределения приведены в 6.4.

Примечание - Настоящий стандарт не ограничивает возможности использования генераторов, отличных от описанных в данном приложении. Однако перед их использованием необходимо убедиться в том, что генерируемые ими последовательности обладают достаточно хорошими статистическими свойствами. Средства тестирования генератора псевдослучайных чисел для равномерного распределения указаны в С.3.2.

С.2 Распределения общего вида

Выборка для любой строго возрастающей одномерной непрерывной функции распределения G_X() может быть получена посредством выборки из равномерного распределения. Для этого:

a) выбирают случайное число из равномерного распределения R(0, 1);

b) определяют , удовлетворяющее условию G_X() = .

Примечание 1 - Требуемая на этапе b) обратная функция  = G_X^-1() может быть найдена аналитически или определена численными методами.

Пример - Входной величине X приписано экспоненциальное распределение с математическим ожиданием х > 0, определяемое плотностью распределения вероятностей g_X() = exp(-/x)/x в области   0 и g_X() = 0 в области  < 0 (см. 6.4.10). Полученная интегрированием плотности распределения вероятностей функция распределения имеет вид G_X() = 1 - exp(-/x) в области   0 и G_X() = 0 в области  < 0. Аналитическое решение дает  = - х ln(1 - ). Этот результат может быть несколько упрощен. Поскольку для случайной переменной Q, подчиняющейся равномерному распределению R(0, 1), случайная переменная (1 - Q) также будет подчиняться равномерному распределению, то можно записать обратную функцию в виде  = - х ln .

Примечание 2 - Численно значение , обычно определяют как точку пересечения с нулем функции G_X() - . Для определения может быть использован алгоритм поиска отрезка, на концах которого функция имеет разные знаки, такой, например, как метод деления отрезка пополам или, что более эффективно, комбинации линейной интерполяции и метода деления отрезка пополам [11].

Примечание 3 - При использовании генератора случайных чисел из равномерного распределения для получения выборки псевдослучайных чисел из другого распределения следует помнить, что выпадение значения  = 0 или  = 1 может привести к сбою генератора. Примером может служить экспоненциальное распределение (см. 6.4.10). Его плотность распределения [формула (9)] для указанных значений не определена. Генератор, описанный в С.3.3, ошибок подобного рода не дает.

С.3 Равномерное распределение

С.3.1 Общие положения

С.3.1.1 Генератор для равномерного распределения является основой для получения псевдослучайных чисел из любого распределения (см. разделы С.2, С.4 и С.6) при наличии соответствующего алгоритма или формулы. При этом качество получаемой выборки из произвольного распределения зависит от качества работы генератора для равномерного распределения и свойств используемого алгоритма преобразования. Таким образом, только генератор, способный воспроизводить выборку из равномерного распределения с хорошими свойствами вместе с хорошим алгоритмом, обеспечивает генерирование псевдослучайных чисел, хорошо согласующихся с заданным распределением.

С.3.1.2 Отсюда вытекает важность тестирования генератора псевдослучайных чисел для равномерного распределения [31]. Если пользователь не уверен в качестве работы генератора, то его не следует использовать до тех пор, пока соответствующее тестирование не будет проведено. В противном случае не исключено получение ошибочных результатов. Рекомендуется использовать средства тестирования согласно [30]. В С.3.3 приведена процедура генерирования псевдослучайных чисел для равномерного распределения, которая успешно прошла указанное тестирование и проста в применении.

С.3.1.3 Параметры процедуры генерирования псевдослучайных чисел, соответствующих равномерному распределению R(0, 1), - входной, выходной, а также являющийся одновременно входным и выходным, - определены в таблице С.1.

Примечание 1 - При задании одного и того же начального числа результатом может быть генерирование одной и той же последовательности случайных чисел. Этот факт является важным элементом регрессионного тестирования, используемого для определения согласованности результатов, полученных программным средством, с результатами предыдущих версий.

Примечание 2 - Некоторые генераторы псевдослучайных чисел при каждом обращении выдают одно случайное, а некоторые - последовательность значений.

Таблица С.1 - Процедура генерирования псевдослучайных чисел для стандартного равномерного распределения

Входной параметр

q - число значений в генерируемой последовательности

Входной/выходной параметр

t - вектор-столбец, элементы которого могут изменяться в ходе выполнения процедуры. Пользователю обычно нет необходимости знать, как изменяются значения этого параметра, но он помогает контролировать выполнение процедуры генерирования случайных чисел. Эти параметры могут быть реализованы как глобальные переменные и не входить явно в формальные параметры процедуры. Один или несколько элементов (определяют начальное число, используемое при следующем обращении к процедуре

Выходной параметр

z - вектор-столбец q значений, составляющих выборку из R(0, 1)

С.3.1.4 Псевдослучайное число х выборки из R(a, b) определяют преобразованием х = а + (b - a)z, где z - элемент выборки из R(0, 1).

С.3.2 Проверка качества случайной выборки

С.3.2.1 Каждый генератор псевдослучайных чисел должен:

a) обладать хорошими статистическими свойствами,

b) предусматривать возможность реализации на любом языке программирования,

c) давать одни и те же результаты для одного и того же начального числа на любом компьютере.

Желательно также, чтобы он был компактным, т.е. простым при реализации. Одним из таких генераторов, близко приближающихся к удовлетворению перечисленных требований, является генератор Вихманна-Хилла [52, 53]. Он использовался во многих приложениях, включая вычисление неопределенности. Однако длина его цикла (количество генерируемых псевдослучайных чисел до их повторения) составляет 2³¹, что сегодня считается для некоторых задач недостаточным. Более того, при его тестировании не по всем критериям были получены положительные результаты [35]. Наконец, этот генератор был разработан для 16-разрядных компьютеров, тогда как сегодня повсеместно используются 32-разрядные и 64-разрядные компьютеры.

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

С.3.2.2 Для комплексной проверки статистических свойств генератор тестируют пакетом программ TestU01 [30]. Этот программный продукт весьма детализирован и включает в себя большое количество статистических тестов, в том числе расширенный пакет краш-тестирования BigCrush. Некоторые генераторы, успешно выдержавшие это тестирование, приведены в списке, составленном Вихманном и Хиллом [54]. В их число входит и усовершенствованный генератор Вихманна-Хилла (см. С.3.3), обладающий следующими свойствами [54]:

a) его просто реализовать на любом языке программирования, он не зависит от побитовых операций, используемых в некоторых генераторах;

b) структура генератора (количество информации, сохраняемой генератором между запросами) невелика и легка в обращении (сравни с параметром t в таблице С.1);

c) он позволяет легко получить несколько последовательностей, необходимых для высокопараллельных приложений, что, вероятно, будет особенностью вычислений неопределенности в будущем;

d) существуют варианты генератора для 32-разрядных и 64-разрядных компьютеров.

С.3.3 Процедура генерирования выборки псевдослучайных чисел из равномерного распределения

С.3.3.1 Как и его предшественник, улучшенный генератор Вихманна-Хилла представляет собой комбинацию конгруэнтных генераторов. Новый генератор сочетает в себе четыре таких генератора, тогда как предыдущая версия сочетала три. Новый генератор имеет период 2¹²¹, приемлемый для любого возможного применения.

С.3.3.2 В таблице С.2 приведено описание улучшенного генератора Вихманна-Хилла для получения псевдослучайных чисел из R(0, 1) для 32-разрядного компьютера.

С.3.3.3 Для 64-разрядных компьютеров шаг а) вычисления, включая (i) и (ii), в таблице С.2 должен быть заменен более простым шагом:

"а) Для j = 1, ..., 4 вычисляют i_j = (a_j x i_j)mod d_j".

С.4 Нормальное распределение

Процедура, описанная в таблице С.3, обеспечивает выбор случайных значений из стандартного нормального распределения N(0, 1) с использованием преобразования Бокса-Мюллера [3]. Случайное значение х из нормального распределения N(, ) вычисляют по формуле х =  + z, где z - случайное значение из N(0, 1).

С.5 Многомерное нормальное распределение

С.5.1 Из всех многомерных распределений наибольший интерес представляет совместное нормальное распределение N(, V), где - вектор математического ожидания размерности n x 1, а V - ковариационная матрица размерности n x n.

С.5.2 Значения случайной переменной из N(, V) [45, 49] могут быть получены путем использования процедуры, описанной в таблице С.4.

Таблица С.2 - Улучшенный генератор Вихманна-Хилла для псевдослучайных чисел из равномерного распределения на интервале (0, 1) для 32-разрядного компьютера

Входной параметр

Нет

Входные/выходные параметры

i1, i2, i3, i4 - целочисленные параметры, которые требуются в качестве входных величин и изменяются в процессе выполнения процедуры. Перед первым обращением к процедуре им присваивают значения от 1 до 2 147 483 647. Между обращениями значения параметров остаются неизменными. Пользователю обычно нет необходимости знать, как изменяются значения этих параметров, используемых в процедуре генерирования псевдослучайных чисел. Данные параметры могут быть реализованы как глобальные переменные и не входить явно в формальные параметры процедуры.

Константы

а, b, с, d - четырехмерные векторы с целочисленными координатами, т.е. а = (а1, а2, а3, а4)T и т.д., такие, что: аT = (11600, 47003, 23000, 33000), bT = (185127, 45688, 93368, 65075), сT = (10379, 10479, 19423, 8123), dT = 2147483123 х (1, 1, 1, 1) + (456, 420, 300, 0). Между обращениями значения констант остаются неизменными.

Выходной параметр

r - псевдослучайное число из R(0, 1).

Алгоритм вычисления

а) Для j = 1, ..., 4: i) Вычисляют ij = aj x (ij mod bj) - cj x [ij/dj], ii) Если ij < 0, то заменяют ij на ij + dj. b) Вычисляют . с) Вычисляют r = w - [w].

Примечание - [w] означает наибольшее целое, не превосходящее w. Запись ij mod bj означает остаток от деления ij на bj.

Таблица С.3 - Генератор псевдослучайных чисел Бокса-Мюллера

Входной параметр

Нет

Выходной параметр

z1, z2 - два случайных значения, полученных независимо из стандартного нормального распределения

Алгоритм вычисления

а) Независимо генерируют случайные числа r1, и r2 из R(0, 1); b) и .

Примечание 1 - Если V - положительно определенная матрица (т.е. все ее собственные значения строго положительны), то множитель Холецкого R единственен [23, страница 204].

Примечание 2 - Если V не является положительно определенной матрицей, то из-за возможных ошибок округления или других причин R может не существовать. Более того, в случаях, когда одно или несколько собственных значений V хотя и положительны, но очень малы, программная реализация алгоритма факторизации Холецкого может оказаться неспособной сформировать матрицу R из-за ошибок округления в арифметике с плавающей запятой. В любой из этих ситуаций рекомендуется вносить в V малые возмущения таким образом, чтобы множитель R для "возмущенной" матрицы V был хорошо определен.

Таблица С.4 - Генератор случайных чисел из многомерного нормального распределения

Входной параметр

n - размерность многомерного нормального распределения - вектор математических ожиданий размерности n x 1 V - ковариационная матрица размерности n x n q - число генерируемых векторов, состоящих из псевдослучайных чисел

Выходной параметр

X - матрица размерности n x q, j-й столбец которой - генерированный случайный вектор из многомерного нормального распределения

Алгоритм вычисления

a) Для матрицы V формируют множитель Холецкого R, т.е. верхнюю треугольную матрицу, удовлетворяющую условию V = RTR (для генерирования q векторов факторизацию матрицы V необходимо выполнить только один раз). b) Генерируют массив Z размерности n x q чисел из стандартного нормального распределения. c) Вычисляют Х =  + RTZ, где 1 - вектор-столбец, состоящий из единиц и имеющий размерность q x 1.

Простая процедура внесения возмущений описана в [49, страница 322] и реализована в генераторе MULTNORM [45].

Примечание 3 - Если V положительно полуопределенная матрица, тогда ее можно представить в виде V = QQ^T, где Q - ортогональная матрица, - диагональная матрица. Тогда матрица ^1/2Q^T может быть использована вместо R для формирования выборки из N(0, V) даже в случае матрицы V неполного ранга.

С.5.3 На рисунке С.1 показано 200 точек, полученных с использованием генератора MULTNORM [45] из N(, V), где

.

Это совместное распределение двух положительно коррелированных величин. Подобные генераторы описаны в других источниках [12].

С.5.4 На рисунке С.1 точки образуют вытянутый наклонный эллипс. Если недиагональные элементы матрицы V заменить нулевыми значениями, то эти точки образуют круг. Если бы элементы главной диагонали были неравны между собой, а недиагональные элементы были равны нулю, то точки образовали бы эллипс с главными осями, параллельными осям графика. Если бы элементы главной диагонали были отрицательными числами (т.е. величины имели отрицательную корреляцию), то главная ось эллипса имела бы не положительный (как на рисунке С.1), а отрицательный наклон.

X - величина 1, Y - величина 2

Рисунок С.1 - Выборка значений из двумерного нормального распределения с положительной корреляцией

С.6 t-распределение

Процедура, описанная в таблице С.5, представляет метод [29], [44, страница 63] генерирования выборки из t-распределения с степенями свободы.

Таблица С.5 - Генератор псевдослучайных чисел для t-распределения

Входной параметр

- число степеней свободы

Выходной параметр

t - выборка из t-распределения с степенями свободы

Алгоритм вычисления

а) Независимо генерируют случайные числа r1, r2 из равномерного распределения R(0, 1). b) Если r1 < 1/2, то вычисляют t = 1/(4r1 - 1) и  = r2/t2; в противном случае вычисляют t = 4r1 - 3 и = r2. c) Если < 1 - /2 или < (1 + t2/)-(+1)/2, то принимают t в качестве выборочного значения из t-распределения; в противном случае повторяют процедуру с шага а).

Примечание - Чтобы стандартное отклонение t-распределения с степенями свободы было конечным, значение должно быть больше двух.