Монте-Карло для чайников: выборы!
Пора возродить рубрику, заброшенную год назад. Тогда мы сообразили на троих и с помощью разноцветных кубиков и крутого Уокера научились очень быстро (за О(1)) выбирать один вариант из N вариантов с произвольно заданными вероятностями.
Сейчас краткий обзор разных методов выбора случайной точки на диске, на окружности, на сфере и на шаре. Во всех этих задачах подразумевается равномерное распределение по поверхности или по объему. К примеру, при выборе точки на единичном диске, вероятность, что она появится на участке с площадью dS должна быть равна dS/π.
[Выбор точки на единичном диске]
Выбор точки на единичном диске
Если ответ нужно получить в декартовых координатах, проще всего применить метод отбраковки:
мы присваиваем
x=rnd*2-1,
y=rnd*2-1,
L=x2+y2,
пока не выполнится условие L≤1
То есть, сначала мы выбрали точку на квадрате, а потом проверили, лежит ли она на диске.
Чуть интереснее задача решается в полярных координатах. Очевидно, что
φ=2*pi*rnd, т.е любое направление равновероятно.
Радиус может принимать значения от 0 до 1, но явно неравномерно. Найдем функцию распределения. При 0≤x≤1 получим
(вероятность, что мы попадем в диск радиусом x, равна отношению площади этого диска к площади всего единичного диска)
Если известна функция распределения случайной величины, мы можем получить ее из стандартного rnd, найдя обратную функцию:
Это правильный ответ, но он требует вычисления квадратного корня, что довольно затратно. Есть любопытный метод обойтись без него. Оказывается, мы получим то же самое распределение, если будем генерить r вот так:
r=max(rnd,rnd)
т.е. брать две случайные величины, равномерно распределенные между 0 и 1, и выбирать бОльшую. Найдем функцию распределения для r. Выражение выше перепишем в виде, принятом в теории вероятностей:
r=max(ξ,η),
ξ,η - две независимые случайные величины, имеющие функцию распределения на [0;1]:
Fξ(x)=Fη(x)=p[ξr(x)=p[r2
(мы пользуемся тем, что вероятности независимых событий перемножаются)
Ч.Т.Д.
[Выбор точки на единичной окружности]
Выбор точки на единичной окружности
В полярных координатах задача элементарная: φ=2π*rnd.
Если нас интересуют декартовы координаты, есть любопытный способ обойтись без математических функций вроде синусов, косинусов или квадратных корней. Начинаем с выбора точки на диске по методу отбраковки (см выше), получаем x0,y0 и L.
А дальше находим:
Это и будет наша точка, лежащая на единичной окружности.
Доказательство.
Сначала найдем точку на окружности более «прямолинейным» методом:
то есть, мы просто взяли произвольную точку на диске и «отнормировали» ее, спроецировали на единичную окружность. Можно написать, что
x1=cosφ
y1=sinφ
для некоторого угла φ. Точка с полярным углом 2φ будет столь же равномерно распределена по окружности, как и с углом φ, найдем декартовы координаты для нее:
x=cos2φ=cos2φ-sin2φ
y=sin2φ=2sinφcosφ
подставляя x1,y1, получаем ответ.
Еще небольшое замечание: никто не запрещает величине L однажды стать равной нулю, тогда получим деление на ноль. Но такое может случиться лишь при работе с весьма замысловатыми генераторами псевдослучайных чисел, которые могут дать 2 раза подряд число строго 0.5! Обычный линейный конгруэнтный генератор на это не способен, у него каждый раз новое число, что, вообще-то говоря, не совсем случайно, иногда (с вероятностью 2-b, где b-разрядность rnd) неплохо бы ему и повториться! Но есть и более хитрые генераторы, где для нахождения текущего числа используется не только предыдущее, но и еще какое-то количество предшествующих чисел, это позволяет получить период повторения последовательности сильно больше 2b. Если использовать их, можно в методе отбраковки на всякий случай проверять не только L≤1, но и L>0.
[Выбор точки на единичном шаре и сфере]
Выбор точки на единичном шаре и сфере
Сюда же можно отнести выбор случайного направления в пространстве, одна из самых частых задач.
Если нужны декартовы координаты, проще всего применить метод отбраковки:
присваиваем
x=rnd*2-1
y=rnd*2-1
z=rnd*2-1
L=x2+y2+z2
до тех пор, пока не выполнится L≤1.
Получим наугад выбранную точку на шаре. Если нужна точка на сфере, то нужно еще добавить:
L=1/sqrt(L)
x=x*L
y=y*L
z=z*L
Этот метод прост для понимания, но снова требует извлечения квадратного корня. Оказывается, и здесь его можно избежать, если применить формализм кватернионов.
Кватернионы сами по себе достойны отдельного «ликбеза», больно интересная штука. Если у комплексного числа две компоненты - действительная и мнимая, то у кватерниона - четыре: одна действительная и три мнимых, при трех разных мнимых единицах i,j,k!
Большая польза кватернионов в том, что они тесно связаны с поворотами в трехмерном пространстве (точно так же, как комплексные числа связаны с поворотами на плоскости). Повороту вокруг оси
(вектор имеет единичную длину) на угол φ можно сопоставить кватернион
И далее можно будет повернуть произвольный вектор a по правилам кватернионного умножения:
Так вот, получить кватернион поворота на произвольное направление и на произвольный угол удивительно легко, это все тот же метод отбраковки, только теперь в 4 измерениях:
присваиваем
a=2*rnd-1,
b=2*rnd-1,
c=2*rnd-1,
d=2*rnd-1,
L=a2+b2+c2+d2,
до тех пор, пока не выполнится L≤1.
Искомый кватернион будет иметь вид
, однако мы не стремимся приводить его к этому виду.
Нам осталось взять какой-нибудь единичный вектор и повернуть его с помощью кватерниона. Вектор должен быть каким-нибудь простым, пусть это будет (1;0;0). Если расписать честно его поворот, получим результат:
Как и обещалось - ни одного квадратного корня, синуса или косинуса! Впрочем, метод весьма и весьма громоздкий, и проблемы начинаются с процедуры отбраковки.
При выборе точки на диске мы с вероятностью π/4≈78% получаем пару (x,y) с первого раза, а среднее количество итераций равно 4/π≈1.27. При выборе точки на шаре вероятность получить точку с первого раза равна (4/3 π)/8=π/6≈52% и среднее количество итераций равно 6/π≈1.91. С 4-мерным шаром дела еще хуже, вероятность с первого раза попасть по нему равна лишь (π^2/4)/16=π^2/64≈31%, среднее количество итераций 3.24, и ведь на каждой нам надо получить 4 случайных числа!
Какой метод быстрее на том или ином «железе» - можно узнать только экспериментально. При генерации непосредственно декартовых координат с применением метода отбраковки используются более простые операции, но может сбиваться конвеер, если с первого раза «ткнуть точку» не получилось.
Сферические координаты.
Проще всего получить «азимут» θ:
θ=2pi*rnd
Если мы похожим способом зададим угол места φ (только от -π/2 до π/2), то получим чрезвычайно частое выпадение полюсов - при φ вблизи π/2 любой азимут приведет нас в одну и ту же небольшую область на сфере.
Чтобы получить правильный метод, найдем плотность распределения fφ(x). Возьмем полосу на сфере, с углом места от φ до φ+dφ. Ее площадь равна 2πcosφdφ,
а полная площадь нашей единичной сферы: 4π. Вероятность, что точка попадет в указанную полосу, равна ½cosφdφ. Плотность распределения, соответственно, равна
fφ(x)=½cosx
Проинтегрировав ее от -π/2 до φ, получим функцию распределения:
Fφ(x)=(1+sinφ)/2
Чтобы из равномерно распределенной на [0;1] случайной величины (rnd) получить φ, нужно просто взять обратную функцию:
φ=arcsin(2*rnd-1)
Если нам нужна точка на сфере, то задача уже решена: мы получили θ и φ - азимут и угол места нашей точки. Для выбора точки на шаре нужен еще радиус. По аналогии с диском:
r=max(rnd,rnd,rnd)
Ясное дело, это далеко не все предложенные методы. Есть методы, основанные на использовании гауссовых случайных величин вместо равномерно распределенных, правда сразу возникает вопрос - откуда взять эти гауссовы случайные числа? И еще есть метод Марсальи для получения случайной точки на сфере - он представляет собой нечто среднее из уже рассмотренных: сначала выбирается случайная точка (x1,y1) на диске методом отбраковки, а потом строится результат:
В итоге нужен всего один квадратный корень, но в отличии от "прямолинейного" метода отбраковки, уменьшается количество итераций (с 1.91 до 1.27) и количество действий на каждой из них, что может сыграть положительную роль.
Литература.
1. Дональд Кнут, Искусство программирования, том 2, глава 3.4.1 - численные распределения. О том, как генерировать случайную величину с заданным распределением, метод отбраковки в общем виде, всякие хитрости вроде взятия максимума или минимума от набора величин.
2. Monte Carlo Ray Tracing - Siggraph 2003 Course 44 Глава 2.3 про выбор точки на диске, сфере, и не обязательно с равномерным распределением
3. Cook, J. M. "Technical Notes and Short Papers: Rational Formulae for the Production of a Spherically Symmetric Probability Distribution." Math. Tables Aids Comput. 11, 81-82, 1957. - это про использование кватернионов для выбора точки на сфере. К сожалению, не смог найти в электронном виде.
4. Marsaglia, G. "Choosing a Point from the Surface of a Sphere." Ann. Math. Stat. 43, 645-646, 1972. - Марсалья рассматривает несколько методов выбора точки на сфере, после чего предлагает свой.
Продолжение следует...
Сейчас краткий обзор разных методов выбора случайной точки на диске, на окружности, на сфере и на шаре. Во всех этих задачах подразумевается равномерное распределение по поверхности или по объему. К примеру, при выборе точки на единичном диске, вероятность, что она появится на участке с площадью dS должна быть равна dS/π.
[Выбор точки на единичном диске]
Выбор точки на единичном диске
Если ответ нужно получить в декартовых координатах, проще всего применить метод отбраковки:
мы присваиваем
x=rnd*2-1,
y=rnd*2-1,
L=x2+y2,
пока не выполнится условие L≤1
То есть, сначала мы выбрали точку на квадрате, а потом проверили, лежит ли она на диске.
Чуть интереснее задача решается в полярных координатах. Очевидно, что
φ=2*pi*rnd, т.е любое направление равновероятно.
Радиус может принимать значения от 0 до 1, но явно неравномерно. Найдем функцию распределения. При 0≤x≤1 получим
(вероятность, что мы попадем в диск радиусом x, равна отношению площади этого диска к площади всего единичного диска)
Если известна функция распределения случайной величины, мы можем получить ее из стандартного rnd, найдя обратную функцию:
Это правильный ответ, но он требует вычисления квадратного корня, что довольно затратно. Есть любопытный метод обойтись без него. Оказывается, мы получим то же самое распределение, если будем генерить r вот так:
r=max(rnd,rnd)
т.е. брать две случайные величины, равномерно распределенные между 0 и 1, и выбирать бОльшую. Найдем функцию распределения для r. Выражение выше перепишем в виде, принятом в теории вероятностей:
r=max(ξ,η),
ξ,η - две независимые случайные величины, имеющие функцию распределения на [0;1]:
Fξ(x)=Fη(x)=p[ξr(x)=p[r2
(мы пользуемся тем, что вероятности независимых событий перемножаются)
Ч.Т.Д.
[Выбор точки на единичной окружности]
Выбор точки на единичной окружности
В полярных координатах задача элементарная: φ=2π*rnd.
Если нас интересуют декартовы координаты, есть любопытный способ обойтись без математических функций вроде синусов, косинусов или квадратных корней. Начинаем с выбора точки на диске по методу отбраковки (см выше), получаем x0,y0 и L.
А дальше находим:
Это и будет наша точка, лежащая на единичной окружности.
Доказательство.
Сначала найдем точку на окружности более «прямолинейным» методом:
то есть, мы просто взяли произвольную точку на диске и «отнормировали» ее, спроецировали на единичную окружность. Можно написать, что
x1=cosφ
y1=sinφ
для некоторого угла φ. Точка с полярным углом 2φ будет столь же равномерно распределена по окружности, как и с углом φ, найдем декартовы координаты для нее:
x=cos2φ=cos2φ-sin2φ
y=sin2φ=2sinφcosφ
подставляя x1,y1, получаем ответ.
Еще небольшое замечание: никто не запрещает величине L однажды стать равной нулю, тогда получим деление на ноль. Но такое может случиться лишь при работе с весьма замысловатыми генераторами псевдослучайных чисел, которые могут дать 2 раза подряд число строго 0.5! Обычный линейный конгруэнтный генератор на это не способен, у него каждый раз новое число, что, вообще-то говоря, не совсем случайно, иногда (с вероятностью 2-b, где b-разрядность rnd) неплохо бы ему и повториться! Но есть и более хитрые генераторы, где для нахождения текущего числа используется не только предыдущее, но и еще какое-то количество предшествующих чисел, это позволяет получить период повторения последовательности сильно больше 2b. Если использовать их, можно в методе отбраковки на всякий случай проверять не только L≤1, но и L>0.
[Выбор точки на единичном шаре и сфере]
Выбор точки на единичном шаре и сфере
Сюда же можно отнести выбор случайного направления в пространстве, одна из самых частых задач.
Если нужны декартовы координаты, проще всего применить метод отбраковки:
присваиваем
x=rnd*2-1
y=rnd*2-1
z=rnd*2-1
L=x2+y2+z2
до тех пор, пока не выполнится L≤1.
Получим наугад выбранную точку на шаре. Если нужна точка на сфере, то нужно еще добавить:
L=1/sqrt(L)
x=x*L
y=y*L
z=z*L
Этот метод прост для понимания, но снова требует извлечения квадратного корня. Оказывается, и здесь его можно избежать, если применить формализм кватернионов.
Кватернионы сами по себе достойны отдельного «ликбеза», больно интересная штука. Если у комплексного числа две компоненты - действительная и мнимая, то у кватерниона - четыре: одна действительная и три мнимых, при трех разных мнимых единицах i,j,k!
Большая польза кватернионов в том, что они тесно связаны с поворотами в трехмерном пространстве (точно так же, как комплексные числа связаны с поворотами на плоскости). Повороту вокруг оси
(вектор имеет единичную длину) на угол φ можно сопоставить кватернион
И далее можно будет повернуть произвольный вектор a по правилам кватернионного умножения:
Так вот, получить кватернион поворота на произвольное направление и на произвольный угол удивительно легко, это все тот же метод отбраковки, только теперь в 4 измерениях:
присваиваем
a=2*rnd-1,
b=2*rnd-1,
c=2*rnd-1,
d=2*rnd-1,
L=a2+b2+c2+d2,
до тех пор, пока не выполнится L≤1.
Искомый кватернион будет иметь вид
, однако мы не стремимся приводить его к этому виду.
Нам осталось взять какой-нибудь единичный вектор и повернуть его с помощью кватерниона. Вектор должен быть каким-нибудь простым, пусть это будет (1;0;0). Если расписать честно его поворот, получим результат:
Как и обещалось - ни одного квадратного корня, синуса или косинуса! Впрочем, метод весьма и весьма громоздкий, и проблемы начинаются с процедуры отбраковки.
При выборе точки на диске мы с вероятностью π/4≈78% получаем пару (x,y) с первого раза, а среднее количество итераций равно 4/π≈1.27. При выборе точки на шаре вероятность получить точку с первого раза равна (4/3 π)/8=π/6≈52% и среднее количество итераций равно 6/π≈1.91. С 4-мерным шаром дела еще хуже, вероятность с первого раза попасть по нему равна лишь (π^2/4)/16=π^2/64≈31%, среднее количество итераций 3.24, и ведь на каждой нам надо получить 4 случайных числа!
Какой метод быстрее на том или ином «железе» - можно узнать только экспериментально. При генерации непосредственно декартовых координат с применением метода отбраковки используются более простые операции, но может сбиваться конвеер, если с первого раза «ткнуть точку» не получилось.
Сферические координаты.
Проще всего получить «азимут» θ:
θ=2pi*rnd
Если мы похожим способом зададим угол места φ (только от -π/2 до π/2), то получим чрезвычайно частое выпадение полюсов - при φ вблизи π/2 любой азимут приведет нас в одну и ту же небольшую область на сфере.
Чтобы получить правильный метод, найдем плотность распределения fφ(x). Возьмем полосу на сфере, с углом места от φ до φ+dφ. Ее площадь равна 2πcosφdφ,
а полная площадь нашей единичной сферы: 4π. Вероятность, что точка попадет в указанную полосу, равна ½cosφdφ. Плотность распределения, соответственно, равна
fφ(x)=½cosx
Проинтегрировав ее от -π/2 до φ, получим функцию распределения:
Fφ(x)=(1+sinφ)/2
Чтобы из равномерно распределенной на [0;1] случайной величины (rnd) получить φ, нужно просто взять обратную функцию:
φ=arcsin(2*rnd-1)
Если нам нужна точка на сфере, то задача уже решена: мы получили θ и φ - азимут и угол места нашей точки. Для выбора точки на шаре нужен еще радиус. По аналогии с диском:
r=max(rnd,rnd,rnd)
Ясное дело, это далеко не все предложенные методы. Есть методы, основанные на использовании гауссовых случайных величин вместо равномерно распределенных, правда сразу возникает вопрос - откуда взять эти гауссовы случайные числа? И еще есть метод Марсальи для получения случайной точки на сфере - он представляет собой нечто среднее из уже рассмотренных: сначала выбирается случайная точка (x1,y1) на диске методом отбраковки, а потом строится результат:
В итоге нужен всего один квадратный корень, но в отличии от "прямолинейного" метода отбраковки, уменьшается количество итераций (с 1.91 до 1.27) и количество действий на каждой из них, что может сыграть положительную роль.
Литература.
1. Дональд Кнут, Искусство программирования, том 2, глава 3.4.1 - численные распределения. О том, как генерировать случайную величину с заданным распределением, метод отбраковки в общем виде, всякие хитрости вроде взятия максимума или минимума от набора величин.
2. Monte Carlo Ray Tracing - Siggraph 2003 Course 44 Глава 2.3 про выбор точки на диске, сфере, и не обязательно с равномерным распределением
3. Cook, J. M. "Technical Notes and Short Papers: Rational Formulae for the Production of a Spherically Symmetric Probability Distribution." Math. Tables Aids Comput. 11, 81-82, 1957. - это про использование кватернионов для выбора точки на сфере. К сожалению, не смог найти в электронном виде.
4. Marsaglia, G. "Choosing a Point from the Surface of a Sphere." Ann. Math. Stat. 43, 645-646, 1972. - Марсалья рассматривает несколько методов выбора точки на сфере, после чего предлагает свой.
Продолжение следует...