Научная статья на тему 'Об оптимальном моделировании дискретных распределений'

Об оптимальном моделировании дискретных распределений Текст научной статьи по специальности «Математика»

CC BY
214
90
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ДИСКРЕТНЫЕ РАСПРЕДЕЛЕНИЯ / МОДЕЛИ ГЕНЕРАТОРОВ ПСЕВДОСЛУЧАЙНЫХ ЧИСЕЛ / СЛОЖНОСТЬ МОДЕЛИРОВАНИЯ / ПРЕДСТАВЛЕНИЕ ЧИСЕЛ В ФОРМАТЕ С ПЛАВАЮЩЕЙ ТОЧКОЙ / DISCRETE DISTRIBUTIONS / MODELS OF PSEUDO-RANDOM GENERATORS / COMPLEXITY OF SIMU- LATION / FLOATING POINT NUMBER REPRESENTATION

Аннотация научной статьи по математике, автор научной работы — Воробьева Н. А., Коробейников А. И., Некруткин В. В.

В статье предполагается, что источником случайности при моделировании дискретных распределений является последовательность независимых равномерно распределенных на множестве {0,...,M − 1} случайных величин. В этом предположении строятся и исследуются оптимальные (в смысле среднего числа обращений к этой модели генератора псевдослучайных чисел) алгоритмы, а также обсуждаются вопросы учета представления чисел в формате с плавающей точкой. Приведены примеры применения полученных результатов к традиционным методам моделирования.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Optimal simulation of discrete distri- butions

An infinite sequence of i.i.d. variables uniformly distributed on the set 0,...,M − 1 is considered as the model of the pseudo-random number generator. Under this assumption, optimal algorithms for simulation of random variables with discrete distributions are constructed and analyzed. The role of floating-point representation of probabilities is discussed. Several examples, related to traditional simulation methods are presented

Текст научной работы на тему «Об оптимальном моделировании дискретных распределений»

УДК 519.245:519.68

Вестник СПбГУ. Сер. 1. 2012. Вып. 3

ОБ ОПТИМАЛЬНОМ МОДЕЛИРОВАНИИ ДИСКРЕТНЫХ РАСПРЕДЕЛЕНИЙ*

H. А. Воробьева1, А. И. Коробейников2, В. В. Некруткин3

I. С.-Петербургский государственный университет, студент, nv.nerely@mail.ru

2. С.-Петербургский государственный университет, канд. физ.-мат. наук, ассистент, anton@korobeynikov.info

3. С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, vnekr@mail.ru

1. Введение и постановка задачи. Пусть X = {xi,..., xi,...} — не более чем счетное множество. Обозначим I = {1,...,N}, если cardX = N < ж, и I = N в противном случае. Рассмотрим задачу моделирования дискретного распределения P, сосредоточенного на множестве X и заданного таблицей распределения

P: ( xi ... xi ... V I е I, (1)

V pi ... pi ... j

предполагая, что источником случайности является последовательность независимых случайных величин ei,..., en,..., каждая из которых имеет распределение

Q : ( 0 ... k ... M - 1 ) , M > 1. (2)

V qo ... qk ... qM-1 J

Традиционно считается, что такая постановка задачи восходит к [1], где обсуждается случай M = N = 2, pi = p2 = 1/2 и qo = qi = q е (0,1), то есть рассматривается моделирование симметричных испытаний Бернулли с помощью несимметричных.

Идеи работы [1] развивались в различных направлениях. Одно из таких направлений, «the Bernoulli factory», изучает свойства функций f : (0,1) ^ (0,1), для которых при помощи испытаний Бернулли с вероятностью успеха q можно промоделировать испытания Бернулли с вероятностью успеха p = f (q), причем число q считается неизвестным. Некоторые современные результаты по этой тематике можно найти в [2-4].

Другое направление — это конструирование и исследование алгоритмов моделирования в достаточной мере произвольных распределений (1) в случае, когда распределение (2) является простым и универсальным. Здесь можно выделить статью [5], где обсуждается оптимальное моделирование распределения P в том случае, когда распределение Q представляет собой симметричное распределение Бернулли. Соответствующие алгоритмы опубликованы в [6] и [7].

Как правило, моделирование распределения (1) с помощью независимых случайных величин ej, имеющих распределение (2), производится методом корневых

* Работа выполнена при частичной финансовой поддержке РФФИ (гранты № 11-01-00769-а и 12-01-00747-а).

© Н.А.Воробьева, А. И. Коробейников, В. В. Некруткин, 2012

диаграмм. Формально нужная корневая диаграмма описывается набором Б цепочек (¿1,...,гг) с ¿т € {0,..., М — 1}, удовлетворяющих условию префикса: если (¿1,. .., V) € Б, то (¿1, .. ., 31,. .., Зи) € Б при всех п > 1 и 31, .. . ,3п € {0,..., М — 1}.

Кроме того, предполагается, что

%=1 (з)

В этом случае мы будем говорить, что корневая диаграмма Б согласована с 'распределением 2, и будем обозначать ее Б = Б(2).

Наконец, пусть корневая диаграмма Б(2) связана с распределением Р при помощи отображения / : Б ^ X такого, что для любого Х£ € X

... ^ = Р (4)

/-1({Х1|)

(конечно, (3) следует из (4)). Тогда будет употребляться термин корневая диаграмма Б(Р, 2, /), порождающая распределение Р.

Положим теперь т = шш{п : (е1,...,е„) € Б}. Тогда, очевидно, (3) означает, что Р(т < то) = 1, равенство (4) переписывается в виде Р(/(£1,... , £т) = Х£) = р^, и общая схема моделирования распределения Р выглядит следующим образом: ищется наименьшее п такое, что (£1,..., £и) € Б, и полагается £ = /(£1,..., £и).

Естественной характеристикой трудоемкости такого метода моделирования является величина Ет — среднее число случайных величин £ ^, необходимое для получения одной реализации случайной величины £. Точное выражение для Ет хорошо известно в теории информации (например, [8, гл. 1 §4]). А именно, если рассмотреть распределение 2б, имеющее носитель Б, и такое, что 2б((¿1... ¿и)) = П^,= 1 , то окажется, что

^(2в) = ЕтЯ2(б), (5)

где Нк обозначает К-энтропию дискретного распределения Как отмечено в [9], равенство (5) непосредственно вытекает из тождества Вальда, примененного к последовательности случайных величин log2(qe ) и моменту остановки т.

Поскольку Н2 (2б) > Н2(Р) (равенство достигается, если / — биекция), а Н2(2) < ^2(М), то из (5) следует неравенство Ет > Н2(Р)/^2(М), носящее в теории информации название обратного неравенства Шеннона для кодирования канала без помех.

Способы построения корневых диаграмм для моделирования распределения Р обсуждаются, например, в [9] и [10]. Однако вопрос об оптимальном (то есть минимизирующем величину Ет) алгоритме, по всей видимости, остается в общем случае открытым.

В настоящей работе мы считаем, что распределение 2 является равномерным на множестве {0,..., М — 1}. Иначе говоря, предполагается, что Цк = 1/М. Для краткости вместо Б(2) и Б(Р, 2, /) мы будем использовать обозначения Бм и Бм(Р, /).

Как уже говорилось, случай М = 2 подробно рассмотрен в [5], где описаны оптимальные корневые Б2(Р, /)-диаграммы и проанализирована соответствующая величина Ет. В принципе, перенос результатов для М = 2 на произвольное М не представляет большого труда, так как логика рассуждений статьи [5] остается прежней.

Тем не менее случай М > 2 представляет особый интерес, связанный с традиционным статистическим моделированием, где источником случайности объявляется последовательность независимых равномерно распределенных на отрезке [0,1] случайных величин (например, [11]), а вероятности р являются произвольными положительными вещественными числами, удовлетворяющими условию ^г р = 1.

На практике же генераторы псевдослучайных чисел (например, линейные конгруэнтные генераторы) порождают последовательность целых чисел, принимающих значения от нуля до М — 1 для некоторого целого М > 1. Поэтому более адекватной моделью источника случайности является последовательность независимых равномерно распределенных на множестве {0,..., М — 1} случайных величин {е^>1, что требует разработки новых алгоритмов моделирования и анализа их трудоемкости.

Чтобы подчеркнуть такую практическую направленность представленных ниже результатов, мы будем интерпретировать реализации случайных величин е^ как результаты обращения к генератору случайных чисел (для краткости — просто к генератору), а величину Ет называть средним числом обращений к генератору, необходимым для моделирования распределения Р.

Кроме того, при машинных вычислениях вероятности р£ имеют представление в виде чисел с плавающей точкой, то есть каждая из них принадлежит конечному подмножеству отрезка (0,1), образующему неправильную решетку.

Этим вопросам и посвящена настоящая статья. Во втором ее параграфе описываются и анализируются оптимальные методы моделирования распределения (1) с произвольными р^. В параграфе 3 рассматриваются соответствующие алгоритмические схемы. В четвертом — учитывается специфика представления чисел с плавающей точкой. Наконец, в последнем параграфе обсуждается применение полученных результатов для анализа традиционных методов моделирования.

2. Оптимальные Бм(Р, /)-диаграммы. Итак, пусть у нас имеется бесконечная последовательность независимых одинаково распределенных случайных величин еi, принимающих с равной вероятностью значения 0,1,..., М — 1 при некотором М > 1. В дальнейшем члены этой последовательности будем называть случайными М-битами.

Рассмотрим некоторую корневую диаграмму Бм (Р, /), порождающую распределение Р, и обозначим ¿г(ш) = еаМ {(¿1,..., ¿т) : /((¿1,..., ¿т)) = х^}. Тогда условие (4) превратится в равенство ^т>1 ¿¿(т)/Мт = р^.

Введем при х € (0, 1) функции

v(M}(x) = Y, {Mmx}/M

г> 0

(здесь {z} обозначает дробную часть числа z) и em(x) = |_MmxJ (mod M), так что x = 2m>i em(x)M-m. Пусть, наконец, £ = f (eb ..., eT).

Предложение 1. 1. Для любой корневой Dm (P, f)-диаграммы и любого m > 0 имеет место неравенство

Р(т > m) > J2 {Mmpi}/Mm. (6)

2. Равенства в (6) достигаются тогда и только тогда, когда ti(m) = em(pi) для всех m и I, при этом величина Ет является минимальной и

Ет = £ v(pi). (7)

i

Доказательство. 1. Очевидно, Р(т < m, £ = xe) = J2i<k<m t^(k)/Mk < что эквивалентно неравенству

]Г mm-kti(k) < LM>ij,

l<k<m

так как сумма в его левой части есть целое число. Отсюда сразу же получаем, что

Р(т > m) = £ Р(т > m, £ = x^ = £ (Р£ - Р(т < m, £ = x^) >

i i

> £ (pi - lmmpij/mm) = £{mmpi}/M

i i

2. Если в (6) достигаются равенства, то

ti(m)+ M £ Mm-k-1ti(k)= LM mpij

1<k<m

для всех I, m. Следовательно, ti(m) = LMmpiJ (mod M) = em(pi). Обратное утверждение очевидно. Представление (7) следует из определения функции v(x) и равенства Er = Em>o Р(т > m)- 1=1

Замечание 1. Результат предложения 1 показывает, что выбор ti(m) = em(pi) обеспечивает не только минимальное среднее число обращений к генератору, но одновременно и минимальные значения вероятностей Р(т > m) при всех m > 1. В частности, при этом достигается максимум вероятности Р(т = 1),

Р(т = 1) = £ LMpiJ/M. (8)

i

Предложение 2. Минимальное значение min Ет удовлетворяет неравенствам

max (1, Hm(P)) < min Ет < hm(P) + M/(M - 1).

Доказательство. Левое неравенство уже обсуждалось во введении. Для доказательства правого неравенства заметим, что

v(x) = £ kek(x)M"k k>1

при x G (0,1), и определим функцию hm(x) равенством hm(x) = —xlogM(x). Ввиду (7) достаточно показать, что v(x) < HM (x) + xM/(M — 1).

Для x G (0,1) найдем такие целые m> 0 и 1 < j < M, что LxM mJ = j. Тогда

Hm(x) + xM/(M - 1) - v(x) = £ ek(x)M"k ( - logM(x) + M/(M - 1) - fc) >

k>1

> jM"т(м/(м-1)-logM(j + 1)) + £ ek(x)M"k(m-k+M/(M-1)-logM(j + 1)) .

m

k>m

(9) 17

Так как М > ] > 1, верно 0 < М/(М — 1)— + 1) < 1, и поэтому первое слагаемое

в правой части (9) положительно, а остальные — отрицательны. Следовательно,

Нм (х) + хМ/(М — 1) — V(х) > ^М-т(М/(М — 1) — ^м + 1)) +

+ ]Т(М — 1)М-*(М/(М — 1) — logм + 1)) — £ (М — 1)М-к(к — т) =

к>т к>т

= М-т (+ 1) (М/(М — 1) — logм + 1)) — М/(М — 1)).

Для окончания доказательства осталось заметить, что у (г/(г — 1) — logz у) > г/(г — 1) при 2 £ (1, +оо) и у £ (1, г]. □

3. Оптимальные алгоритмы моделирования. В этом разделе мы воспользуемся результатами предложения 1 для построения оптимальных (M, P)-алгоритмов, преобразующих случайные M-биты в случайную величину £ с распределением (1), где

pi =$3 em(pi)M

m> 1

а em(x) = |MmxJ (mod M).

Как уже говорилось, при M = 2 такие алгоритмы рассматривались в [6] (при конечном N) ив [7]. В принципе, идеи [6] и [7] несложным образом переносятся и на случай M > 2. Поэтому мы приводим соответствующие алгоритмы без специальных объяснений.

Предполагается, что на вход алгоритмов поступают целые числа em(pi). Вызов функции Um () означает «получить реализацию случайной величины, имеющей равномерное распределение на множестве {0,..., M — 1}». Результат моделирования обозначен через £.

Алгоритм 1: Общий, cardX = N < то.

2 for m > 1 do

3 e ^ Um(); j ^ M x j — e ;

4 for I ^ 1 to N do

5 j ^ j — em (pi) ;

6 if j < 0 then (£ ^ x£; STOP)

7 end

8 end

Конечно, в некоторых частных случаях алгоритм 1 существенно упрощается. Например, если Р представляет собой распределение Бернулли с параметром р, то алгоритм 1 можно записать в следующем виде (см. [12] при М = 2).

Алгоритм 2: Распределение Бернулли.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

1 m ^ 0 ;

2 repeat

3 m ^ m + 1; е ^ Um() ;

4 if е < em(p) then (£ ^ 1; STOP)

5 until е > em(p) + em(1 - p) ;

6 e ^ о

При больших N алгоритм 1 весьма трудоемок. Кроме того, для его реализации требуется умение работать с очень большими целыми числами.

Действительно, пусть N = Мк + 1 для к > 0. Рассмотрим такой набор вероятностей рх,... , что вк+1 (р) = 1 при г = 1,..., N — 1 и ) = М — 1. Несложно видеть, что ¿=1 Рг = 1. Теперь при работе алгоритма 1 максимальное значение переменной ] в строке 3 составит Мк и будет достигаться на последовательности £ = 0,..., 0 (к раз). Ситуация, аналогичная представленной, может возникать в том случае, когда набор вероятностей рх,... содержит как достаточно большие, так и достаточно маленькие значения одновременно.

Следующий алгоритм (см. [7] для случая М = 2) основан на методе уточнения, предложенном в [5]. Он использует алгоритм 2 и формально может быть использован и при I = N. Положим

91 = £р = £ ет(®)М-т. (10)

г>1 т>1

Алгоритм 3: Общий, метод уточнения.

1 m ^ 0; I ^ 1 ;

2 while m > 0 do

3 repeat

4 m ^ m + 1; e ^ Um() ;

5 if e < em(pi) then (£ ^ xi; STOP)

6 until e > em(pi) + em(qi) ;

7 repeat

8 e ^ e - em(pi); I ^ I +1 ;

9 if e < em(pi) then (£ ^ xi; STOP)

10 until e < em(pi) + em(qi) ;

11 end

4. Учет представления чисел с плавающей точкой. В предыдущих разделах предполагалось, что вероятности pi являются вещественными числами с известными коэффициентами en(pi) разложения по степеням M-n. Если это не так, то при компьютерном моделировании числа pi получают представление с плавающей точкой, которое определяется двумя целыми положительными параметрами B и S — количеством двоичных бит соответственно в порядке и мантиссе.

Всюду далее предполагается, что число с плавающей точкой имеет представление в соответствии с общеупотребительным стандартом IEEE 754-2008 [13], при этом для простоты изложения игнорируются денормализованные числа. Заметим, что их учет в представленных ниже алгоритмах не представляет сложности.

Тогда все числа, принадлежащие промежутку (0,1), имеют вид

хд-,к = 2— (1 + к 2-й) , (11)

где к = 0,..., 25 — 1, ] = 1,..., Ь и Ь = 22 -2. Таким образом, минимальное из чисел х^,к равно 2-Ь, максимальное равно 1 — 2-5, а всего таких чисел Ь25. Отметим, что в случае одинарной точности В = 8 и Б = 23, для двойной точности В = 11 и Б = 52, а для чисел с расширенной точностью В = 15 и Б = 64. Для наиболее часто употребляемых чисел с двойной точностью минимальное из чисел х^к составляет 2-1022 « 2.2 х 10-308 (при использовании денормализованных чисел — 4.9 х 10-324), чего обычно хватает для большинства приложений.

Представление вероятностей р в виде (11) вносит свою специфику в реализацию алгоритмов 1-3. Обсудим этот вопрос в случае, когда М = 2^. Сначала получим явное выражение для М-битов е„(х^,к) чисел х^к. Положим п* = Г.?/^] и п* = + Б)/Щ]. Ясно, что п* < |~2В/Щ и п* < |"(2В + Б)/Щ. Тогда из формулы (11) непосредственно следует, что е„(х^к) = 0 при п<п* и п>п*,ав остальных случаях

(+ [k/2S+j-"*dJ при n : e„(xj,fc) = \ [k/2S+j_ndJ (mod 2d) при n* < n < n*, (12)

(k2n*d_S_j (mod 2d) при n = n*

Отсюда сразу видно, что затраты памяти на хранение таблицы значений en(xj,fc) составляют n* — n* + 1 < [S/dJ + 2 d-битовых слов. Тем самым в худшем случае мы здесь тратим на 2 d-битовых слова больше, чем необходимо для хранения одного числа Xj fc в представлении с плавающей точкой.

Вернемся к алгоритму 1 и предположим, что вещественные числа pi получают представление с плавающей точкой (напомним, что мы пользуемся вариантом стандарта IEEE 754-2008 [13] без учета денормализованных чисел).

Если все pi уже принадлежат решетке (11), то алгоритм 1 остается точным. Если же это условие не выполнено, то числа pi подвергаются автоматическому округлению и превращаются в числа p*, имеющие вид (11). В случае, когда среди вероятностей pi есть очень маленькие числа, результатом их округления может оказаться 0, если есть очень близкие к 1, то они могут округлиться до 1. Будем считать, что такие эффекты отсутствуют, то есть min pi > 2_L-1 и max pi < 1 — 2_S_1.

Кроме того, после округления условие ¿=i pi = 1, вообще говоря, нарушается, и алгоритм 1 перестает быть корректным. Чтобы избежать подобных неприятностей, можно получить коэффициенты (12) для всех p*, кроме последнего pV, обозначить pv = 1 — X^iLi1 p* и вычислить коэффициенты en(pv) через em(p*), 1 < I < N. Тогда алгоритм 1 будет оптимальным и точным для вероятностей p*,... , pV _i,pv.

Это реализуемая, но громоздкая операция. Покажем, как она решается для распределения Бернулли.

Если p = имеет представление (11), то из формулы (12) несложно получить, что en(1 — p) = 0 при n > n* и

2d — en(xj,fc) — 1 при n < n 2d — en(xj,fc) при n = n

*

en(1 — )H od Г _ V (13)

Алгоритм 4: Распределение Бернулли, плавающая точка.

1 for n =1 to n* do

2 e ^ Um() ;

3 if e < e„(p) then (£ ^ 1; STOP) ;

4 else if (e = 2d - 1 or n = n*) then (£ ^ 0; STOP)

5 end

6 e ^ 0

Теперь, замечая, что

(2d - 1 при n < n*, 2d при n = n*,

0 при n > n*,

перепишем алгоритм 2 в виде, учитывающем представление числа p с плавающей точкой, предполагая, что числа en(p) уже заранее сосчитаны согласно формуле (12). Результат представлен в виде алгоритма 4.

В случае алгоритма 3 ситуация становится более сложной, так как для его корректности требуется выполнение N равенств вида l=1 pj = qi, l = 1,..., N. Задача, однако, значительно упрощается, если задавать распределение P не вероятностями pi, а суммами qi с 1 < l < N — 1.

При этом, конечно, нужно позаботиться, чтобы числа qi не «слипались» (и не превращались в 0 или 1) при их округлении до чисел q*, принадлежащих решетке (11). Для этого достаточно выполнения неравенства min pi > 2-S-1.

Пусть теперь числа qi принадлежат решетке (11) и коэффициенты en(qi) уже вычислены. Для реализации алгоритма 3 необходимо еще вычислить коэффициенты en(pi). При l = N это просто, так как pw = qw-1. При l =1 можно воспользоваться формулой (13). Для остальных l, как нетрудно видеть, вычисления можно производить согласно следующему алгоритму.

Алгоритм 5: Вычисление коэффициентов en(pi), 2 < I < N — 1.

1 for I = 2 to N - 1 do

2 n ^ n*(qi) + 1; c„ ^ 0;

3 while n > 1 do

4 n ^ n — 1; T ^ en(qi_i) — e„(qi) — c„+i;

5 if T > 0 then (e„(pi) ^ T; c„ ^ 0) else (e„(pi) ^ 2d + T; c„ ^ 1)

6 end

7 end

Заметим, что алгоритм 5 будет охватывать и крайние случаи I =1, I = N, если положить е„(до) = е„(д№) = 0.

Оценим теперь дополнительную память, необходимую для хранения коэффициентов еп(р1). Обозначим через Е порядки чисел 91. Заметим, что из формулы (10) следует Е < Е1+1 и все Е < 0. Так как р = 91-1 — 91, для точного представления р достаточно мантиссы длиной $ + Е1— — Е + 1 двоичных бит. Принимая во внимание,

что Eo = 0 и pn = -1, получаем, что общие затраты памяти на вычисление всех значений pi, I = 1,..., N, составляют не более N(S + 1) + En-i бит.

Таким образом, алгоритм 3 (в отличие от алгоритма 1) оказывается пригодным к реализации в арифметике с плавающей точкой, хотя для этого могут потребоваться большие затраты на предварительные вычисления и дополнительную память.

5. Обсуждение. Многие традиционные методы моделирования дискретных случайных величин используют преобразование вида £ = f (а), где а — случайная величина, равномерно распределенная на [0,1], а f — специальным образом подобранная функция. Например, для моделирования равномерного распределения на множестве {1,...,N} используют моделирующую формулу £ = [Na], а для геометрического распределения с параметром p £ (0,1) —преобразование £ = |_ln(a)/ln(1 — p)J .

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

Попробуем описать ситуации, когда эти погрешности могут оказаться существенными. Предполагая, что источником случайности является последовательность случайных M-битов, мы можем сосчитать вероятность того, что для оптимального моделирования распределения (1) потребуется ровно одно обращение к генератору, см. формулу (8). Если эта вероятность окажется существенно меньше 1, то любой метод моделирования распределения (1), использующий ровно одно обращение к генератору случайных (и, следовательно, псевдослучайных) чисел, будет неточным.

Приведем примеры, относящиеся к равномерному на множестве {1,..., N} и геометрическому распределениям. В первом случае Р(т = 1) = Л|Л—1J с Л = N/M, что равно нулю при N > M. При N < M ситуация другая. Например, при типичном выборе M = 232 для N = 2k < 232 требуется ровно одно обращение к генератору случайных чисел, в то время как при N = 2k + j с 1 < j < 2k

Р(т = 1) = 1 — 2—32+k + j2—k — j2—32,

что дает Р(т = 1) ~ 1 — 2—32+k при маленьких j и больших k. Тем самым случай N, близких к M, требует внимательности и осторожности.

Вероятности одного обращения к генератору. Геометрическое распределение

d к

16 9 10 11 12 13 14 15

32 25 26 27 28 29 30 31

64 57 58 59 60 61 62 63

Р(т = 1) 0.97 0.95 0.92 0.86 0.75 0.59 0.35

Для геометрического распределения удобное для теоретического анализа представление правой части (8), вообще говоря, отсутствует. В таблице приведены результаты компьютерных вычислений вероятности Р(т = 1) для М = 2^ и различных р вида 2-к.

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

Приближенное равенство вероятностей Р(т = 1) в таблице при разных d, k с одинаковой разностью d—k объясняется тем, что при маленьких p число тех значений j, при которых [Mp(1 — p)j J = n > 0, примерно равно ln(1 + 1/n)/p.

Литература

1. von Neumann J. Various techniques used in connection with random digits // Nat. Bur. Stand. Appl. Math. Ser., 1951. Vol.12. P. 36-38.

2. Nacu S., Peres Y. Fast simulation of new coins from old // Annals of Applied Probability, 2005. Vol. 15, 1A. P. 93-115.

3. Holtz O., Nazarov F., and Peres Y. New coins from old, smoothly // Constructive Approximation, 2011. Vol.33, N3. P. 331-363.

4. Latuszynski K. и др. Simulating events of unknown probabilities via reverse time martingales // Random Structures and Algorithms, 2011. Vol.38, N4. P. 441-452.

5. Кнут Д., Яо Э. Сложность моделирования неравномерных распределений // Кибе-ренетический сборник, новая серия. Вып. 19. М.: Мир, 1983. С. 97-158.

6. Походзей Б. Б. Преобразование случайных битов в случайные величины с конечными дискретными распределениями // Вестн. Ленингр. ун-та, 1983. Сер. 1, №13. С. 31-36.

7. Походзей Б. Б. Преобразование случайных битов в случайные величины с произвольными дискретными распределениями // Вестн. Ленингр. ун-та, 1985. Сер. 1, №1. С. 39-43.

8. Чисар И., Кёрнер Я. Теория информации. Теоремы кодирования для дискретных систем без памяти. М.: Мир, 1985. 395 с.

9. Romik D. Sharp entropy bounds for discrete statistical simulation // Statistics & Probability Letters, 1999. Vol. 42. P. 219-227.

10. Han T. S., Hoshi M. Interval algorithm for random number generation // IEEE Transactions on Information Theory, 1997. Vol.43, N2. P. 599-611.

11. Devroye L. Non-Uniform Random Variate Generation, Springer-Verlag, New York; Berlin; Heidelberg; Tokyo, 1986. 817 p.

12. Походзей Б. Б. Оптимальный метод моделирования распределения Бернулли // Изв. АН СССР. Техн. киберн. 1982. Вып. 1. С. 207-208.

13. IEEE Standard for Floating-Point Arithmetic // IEEE Std 754-2008. P. 1-58.

Статья поступила в редакцию 26 апреля 2012 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.