Из равномерного в Гаусса

[Alhimik]

Всем прива!!!! Дело в следующем. Есть массив случайных чисел( $X= \{x_0, x_1, x_2, ..., x_{N-1} \}$ , где N - четное ) распределенных по равномерному закону на интервале $x \in [0;1] $. Есть мнение, что из этого можно получить массив случайных чисел, распределенных по закону Gauss'a с помощью следующего перехода:
$y_j= \sigma \sqrt{-2 \ln x_{2i}} \cos {(2 \pi x_{2i+1})} $
при этом последний массив в два раза меньше, $j=0...\frac N 2 -1?, i=0..N-1$.

Насколько эта информация отвечает действительности? Если да - то как это доказать? Математическое ожидание в новом массиве Y будет таким же, как и в случае равномерного (m=0.5)? Если да, то как его менять.

И еще, может кто-то знает, почему в С++ нет встроенного генератора случайных чисел по закону Gauss'a, а в Delphi - есть?

[venco]

Alhimik в сообщении #381817 писал(а):

$y_j= \sigma \sqrt{-2 \ln x_{2i}} \cos {(2 \pi x_{2i+1})} $
при этом последний массив в два раза меньше, $j=0...\frac N 2 -1?, i=0..N-1$.

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

Alhimik в сообщении #381817 писал(а):

Насколько эта информация отвечает действительности? Если да - то как это доказать?

Соответствует.
Проинтегрировать.

Alhimik в сообщении #381817 писал(а):

Математическое ожидание в новом массиве Y будет таким же, как и в случае равномерного (m=0.5)?

Нет. Матожидание будет 0. Если нужно другое, добавьте константу.

Кстати этот метод не очень хорош. Есть более быстрый и численно устойчивый генератор. См., например, здесь: http://www.taygeta.com/random/gaussian.html

[ewert]

Alhimik в сообщении #381817 писал(а):

Если да - то как это доказать?

Идея состоит в том, что если пары случайных величин $(x_i,y_i)$ распределены независимо и стандартно-нормально (прошу прощения за вульгарность, не величин, конечно, а их реализаций, но лень наводить порядок), то после перехода к полярным координатам полярные углы распределены равномерно, а радиусы -- по показательному закону (последнее, видимо, и имелось в виду под "проинтегрировать"). И наоборот. Показательное же распределение (в отличие от нормального) моделируется равномерным уже элементарно. Отсюда и алгоритм: мы фактически моделируем не одну нормально распределённую случайную величину, а пары независимых. Вы же выковыряли из этих пар только иксовые компоненты, на что, конечно, имели право, но и игрековые тоже вполне можно было оставить (это про намёк насчёт синусов).

[Евгений Машеров]

Alhimik

1. Соответствует. Доказательство есть у Кнута, Ермакова и т.п. В общем, любая книга, где как-то поминается Монте-Карло.
2. Из двух равномерных получается пара нормальных, вторая - по формуле с синусом вместо косинуса. При этом они, несмотря на использование одинаковых равномерных, независимы (для параноиков - набираем массив и перед использованием, употребляя третий генератор, тасуем; для суперпараноиков - четвёртым генератором, скажем, на сдвиговых регистрах, ещё и знаки меняем).
3. Матожидание будет нулевым, дисперсия единичной. Для желаемых параметров надо сперва умножить на требуемое СКО, а затем прибавить "матожидание".
4. Почему нет встроенного? Это науке неизвестно. Может, оттого, что для С много доступных библиотек, и считают, что употребят их, и не надо встраивать.

venco

Существенно лучше использующего тригонометрию алгоритм с отбрасыванием был лет 30 тому. Когда синус считался посредством нетривиальной подпрограммы, а умножение было аппаратное. Что я ещё застал. Как появилась аппаратная тригонометрия (IEEE), преимущества как-то скукожились. А насчёт численной устойчивости... Не вполне уверен.

[Хорхе]

Этот и множество других методов есть в книге Devroye Non-Uniform Random Variate Generation.

[Евгений Машеров]

Вот ещё:
George Marsaglia. Random Number Generators
http://dsp-book.narod.ru/randgen.pdf

[Alhimik]

Евгений Маршалов, а как именно называются эти источники с доказательствами. Google выдаёт на этих товарищей: Кнут "Искусство программирования", а Ермаков "Курс высшей математики для экнономистов". Это оно?

Ewert, не понял причем тут полярная система координат. Если $x_i $ и $y_i$ независимые случайные числа по равномерному закону на интервале [0;1]. В таком случае радиус и угол соответственно: $r_i= \sqrt{x_i^2+y_i^2}$, $\alpha_i = \arctg {\frac {y_i} {x_i}} $. Я правильно представляю? И что именно тут интегрировать?

[ewert]

Пусть $(X,Y)$ -- пара независимых стандартных нормальных величин, т.е. их совместная плотность вероятности $f_{XY}(x,y)=\mathrm{const}\cdot e^{-(x^2+y^2)/2}$ (нормировочные константы не выписываем, т.к. они принципиально не важны). Тогда $X=R\cdot\cos\Phi$ и $Y=R\cdot\sin\Phi$, где, очевидно, случайные величины $R=\sqrt{X^2+Y^2}$ и $\Phi$ независимы, причём $\Phi$ распределена равномерно на отрезке $[0;\,2\pi]$ (это можно и формально обосновать, но стоит ли возиться). Вероятность же попадания $R$ в промежуток $[r;\,r+dr]$ получается интегрированием совместной плотности по бесконечно узкому кольцу радиуса $r$ и ширины $dr$, т.е. равна $\mathrm{const}\cdot e^{-r^2/2}\cdot2\pi r\,dr=2\pi\cdot\mathrm{const}\cdot e^{-r^2/2}d(r^2/2)$, т.е. случайная величина $T={1\over2}R^2$ распределена по стандартному показательному закону с плотностью вероятности $f_T(t)=e^{-t}\ (t>0)$. Тогда $T=-\ln U_1$ и, значит, $R=\sqrt{-2\,\ln U_1}$, где $U_1$ -- равномерно распределённая на $[0;\,1]$ величина. Ну а $\Phi=2\pi\,U_2$, где $U_2$ -- тоже стандартная равномерная величина, независимая от $U_1$. Вот и получается, что

$X=\sqrt{-2\,\ln U_1}\cdot\cos(2\pi\,U_2),$
$Y=\sqrt{-2\,\ln U_1}\cdot\sin(2\pi\,U_2).$

[Евгений Машеров]

А. По ссылкам:
1. Кнут Дональд. Искусство программирования для ЭВМ. Т.2, "Получисленные алгоритмы", гл. 3 (точнее не скажу, надо смотреть по конкретному изданию, но доказательство там есть)
2. Ермаков С.М. Метод Монте-Карло и смежные вопросы (тж. Метод Монте-Карло в вычислительной математике или Избранные алгоритмы метода Монте-Карло, в соавт.). Точные ссылки на страницу не могу дать, тем более много изданий).
3. Соболь И.М. Численные методы Монте-Карло. М., Наука, 1973. Раздел 3.2.2., стр. 63.
Б. По методу:
Сначала делается разумный вывод, что случайных чисел с N(0,1) нам понадобится много, и генерировать сразу два оправдано. Затем очевидно, что два независимых стандартно-нормальных с.ч. это то же, что одно двумерное нормальное с нулевой корреляцией. Если нарисуем плотность его распределения, видим, что линии равной плотности - окружности, так что естественно перейти к полярным координатам (при том, что декартовы - это наши искомые числа), аргумент будет равномерно распределён на окружности, остаётся сгенерировать модуль, и после приведенных выше выкладок обнаруживается, что модуль имеет показательное распределение, а для него обратная функция распределения не простая, а очень простая (с.ч. с любым законом распределения можно сгенерировать, если равномерное с.ч. подставить в функцию, обратную к требуемой функции распределения - но обычно такая функция очень сложна; а тут логарифм и всё). Поэтому генерирует два случайных числа - одно даёт модуль, другое аргумент, затем переводим в декартовы координаты - получаем два нормальных с.ч. Если синус и косинус вычислять дорого (сейчас не очень, а когда-то было весьма накладно; или сейчас, но на каких-нибудь сигнальных процессорах, где тригонометрические функции аппаратно не поддерживаются), то есть вариант, когда, сгенерировав два числа R(-1;1) вначале, смотрят, не превышает ли сумма их квадратов w единицу, и если так - отбрасывают и повторяют генерацию (то есть ищут пару случайных чисел таких, что заданная ими точка внутри единичной окружности), а если не превышает - вычисляют два нормальных (при этом логарифмируют эту сумму квадратов w, а вместо синуса и косинуса употребляют $R_1/w$ и $R_2/w$ соответственно. Требуемое число равномерных чисел тут вырастает в $4/\pi($=1.273 раза, зато нет тригонометрии.
В. Это лишь один из приёмов генерации. Обзоры есть в указанных книгах. И вообще, господа, истязайте себя Кнутом!
Г. Спасибо за повышение в чине меня (или моего папы?), но я Машеров.

[ewert]

Евгений Машеров в сообщении #382396 писал(а):

а тут логарифм и всё

Не совсем всё: надо ещё учесть, что величины $U$ и $1-U$ (где $U$ стандартно-равномерно) распределены одинаково.

Евгений Машеров в сообщении #382396 писал(а):

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

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

[Евгений Машеров]

Я имею в виду, что для вычислений нужно лишь логарифм вычислять, для доказательства надо действительно это учитывать.

[Alhimik]

Да, у Кнута толковые томики и внятным языком, а книженция Ермакова по стилю оформления мне сразу напомнила стиль Берманта, Арамановича "Краткий курс математического анализа" - моего самого первого толкового учебника по высшей математике.
В общем, спасибо всем за ответы.

[erwins]

в c++0x есть

[dim1r]

Alhimik в сообщении #381817 писал(а):

$y_j= \sigma \sqrt{-2 \ln x_{2i}} \cos {(2 \pi x_{2i+1})} $

странно, а если $x_{2i} \rightarrow 0 $ , то $ y_j \rightarrow \infty $ ?

[--mS--]

Разумеется. А что же в этом странного?