nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Генерация гауссовых случайных величин и "наглядное пособие" по отбраковке

В прошлом году я преподавал студентам магистратуры "Алмаза" предмет своего научрука "Корреляционные методы в пассивной локации и ориентации" и счел нужным буквально 2-3 пары посвятить методам Монте-Карло, потому что студенты Физтеха и Бауманки, которые попадали туда, имели пробел в этой области. Они хорошо программировали, и у них был предмет "вычислительная математика", где давались базовые знания по численному решению нелинейных уравнений, систем линейных уравнений и уравнений в частных производных, но при этом ни слова не говорилось о Монте-Карло. Теория вероятностей была на первом курсе, но к 5-му успешно забывалась.

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

Большой проблемой было время - надо было "переть напролом", не сильно залипая на частностях.

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


На каждой итерации мы генерируем случайный вектор, находящийся в определенной области, чаще всего в квадрате, кубе или гиперкубе. С вероятностью p он попадает внутрь "полезной", интересующей нас области, и тогда мы прерываем цикл. Случайную величину, характеризующую число итераций, назовем ξ.

Получается, что с вероятностью p произойдет 1 итерация, с вероятностью (1-p)p - две итерации, с вероятностью (1-p)2p - три итерации, и т.д, получается, f(n)=P[ξ=n]=(1-p)n-1p.
Математическое ожидание равно
eq15.png
Не ахти какая сложная задачка. Заменяем (1-p)=q, выносим p за знак суммы:
eq16.png
Вспоминаем сумму бесконечной геометрической прогрессии:
eq17.png
и замечаем, что нужная нам сумма ряда является производной по q суммы бесконечной геометрической прогрессии:
eq18.png
Ответ совсем прост: 1/p.

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

Но все это из пушки по воробьям, данная задача должна решаться в 1 строчку и быть понятной "на ощупь". Для этого надо вспомнить (или познакомиться) с еще одним определением математического ожидания для положительных целочисленных случайных величин:
eq19.png
Можно это равенство доказать, выразив P[ξ≥n] через сумму P[ξ=n] от n до ∞ и перегруппировав суммы, а можно посмотреть на картинку:
exp.png
У нас есть отрезок [0;1], на котором друг за другом отложены вероятности того, что случайная величина равна 1,2,3 и т.д. Площадь изображенной фигуры равна математическому ожиданию случайной величины. Как видно, посчитать ее можно разными способами, разбивая на горизонтальные или на вертикальные полоски, они-то и приведены в формуле выше.

Так вот, для решения нашей задачки надо определенно использовать нижний способ. Вероятность, что количество итераций будет равно или больше n, равна (1-p)n-1. При суммировании от 1 до ∞ получим

E[ξ]=1+(1-p)+(1-p)2+(1-p)3+... = 1/(1-(1-p)) = 1/p.

Все! А освободившееся время используем, чтобы получше разобраться в следующем методе.

Генерация Гауссовых случайных чисел с нулевым средним и единичным среднеквадратичным
Может быть, это самая распространенная задача в методах Монте-Карло, ведь Гаусс - куда более естественное распределение, чем равномерное на [0..1]. Методов придумано великое множество (см. Кнута 2-й том), один другого веселее. Я в свое время попробовал реализовать два из них - метод полярных координат и метод прямоугольника-клина-хвоста. Последний - сущее мучение для реализации - сначала надо подготовить прорву массивов, да и основной алгоритм весьма и весьма объемный. Что самое обидное - никакого выигрыша по сравнению с методом полярных координат не получилось. Наверное, он был на старых компьютерах, а сейчас гораздо эффективнее "перемалывать числа", пусть даже с логарифмами и квадратными корнями, чем проходить ветвления и обращаться к большим таблицам с риском промахнуться мимо кэша, а время от времени все равно попадая на трансцендентные функции.

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

присваиваем
x=rnd*2-1
y=rnd*2-1
S=x2+y2,
пока не выполнится S≤1
(и снова мы неявно подразумеваем, что x и y при псевдослучайном rnd не могут одновременно равняться нулю и такое условие можно не проверять)

Теперь можно одним махом получить две независимые нормально распределенные случайные величины:
eq20.png

Докажем, что это так. Представим x,y в полярных координатах:

x=Rcosφ,
y=Rsinφ

φ - равномерно распределенная на [0;2π) случайная величина, а R имеет функцию распределения на [0;1):

Fr(x)=x2 (см. "Монте-Карло для чайников: выборы!"), кроме того, R и φ независимы.

Выразим через φ и R значение S и результирующие величины V1, V2:
eq30.png
Величина S является равномерно распределенной на [0;1), ведь
eq23.png
(а еще лучше вспомнить, как мы получали распределение для R, через площадь)

Величина φ' равномерно распределена на [0;2π], осталось найти распределение для R':
eq31.png
Уже выглядывает что-то похожее на Гаусса! Осталось немножко. Записываем плотность распределения для R' и φ':

eq27.png

Перемножив их, мы получим такую величину:
eq26.png
т.е мы записали вероятность попадания в сегмент R' < r < R'+dr, φ' < φ < φ'+dφ. Справа в скобке получаем площадь этого сегмента, ее же можно выразить и в декартовых координатах как dxdy, это позволяет записать двумерную плотность распределения для искомых величин V1, V2:

eq29.png,
а это ни что иное, как нормальное распределение для двух некоррелированных случайных величин с нулевым средним и единичным среднеквадратичным отклонением.

Если мы умеем генерировать такие нормальные случайные величины V, не составит труда получить нормальную случайную величину со средним a и среднеквадратичным отклонением σ:

W=a+σV.
Tags: Монте-Карло для чайников, математика, моделирование
Subscribe

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Потёмкинская деревня - 2

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

  • Очередная несуразность в единицах измерения

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

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 6 comments