Комментарии 41
Подход с треугольниками неаккуратный- он дает распределение, очень близкое к случайному, но отличающееся от него в самой близи от границы круга. дело в том, что Ваш узкий треугольник- это все-таки треугольник, и между его основанием и окружностью есть небольшой участок, в котором точка уже за треугольником, но еще внутри круга. Как следствие- такие точки чуть-чуть отличаются по вероятности возникновения от "нормальных". На практике при малых тетах и небольших выборках это невозможно отследить, но для моделирования чего-нибудь методом Монте-Карло- вполне может дать заметные чудеса.
Выполнив каждую функцию 3 142 592 раза
Мне кажется, там должно быть 3 141 592 :) Интересная статья, спасибо!
Еще бы замеры скорости привести, потому что у меня есть стойкое ощущение, что для быстрого генератора случайных чисел (как правило, это вихрь Мерсенна) будет быстрее отбросить "неправильную" пару случайных чисел, чем рассчитывать значения синуса и косинуса.
Дак они же приведены в самом конце. Правда, судя по всему, на чистом питоне, что несколько снижает ценность такого замера.
Да, в Python для подобных задач надо использовать numpy. Тогда и миллионы точек будут вычисляться за миллисекунды.
def sqrt_dist(n):
theta = np.random.rand(n) * 2 * np.pi
r = np.sqrt(np.random.rand(n))
return r * np.cos(theta), r * np.sin(theta)
%timeit x, y = sqrt_dist(3141592)
139 ms ± 9.76 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Можно ускорить еще почти в два раза, просто поставив декоратор numba.njit:
import numpy as np
from numba import njit
def sqrt_dist(n):
theta = np.random.rand(n) * 2 * np.pi
r = np.sqrt(np.random.rand(n))
return r * np.cos(theta), r * np.sin(theta)
@njit
def sqrt_dist_nb(n):
theta = np.random.rand(n) * 2 * np.pi
r = np.sqrt(np.random.rand(n))
return r * np.cos(theta), r * np.sin(theta)
%timeit x, y = sqrt_dist(3141592)
170 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit x, y = sqrt_dist_nb(3141592)
91.4 ms ± 715 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Только участок BD в треугольнике на самом деле является частью окружности.
Т.е. в данном алгоритме существуют области (предельно малые - но они есть), где точек принципиально не будет.
Я правильно понимаю?
мне тоже так померещилось сначала. но выбрав треугольник безконечно узким мы эту проблему в пределе полностью устраним.
IMHO - численными методами проблему полностью не устраним. "В пределе" - это математика.
А тут уменьшая "толщину" треугольников мы увеличиваем покрытие круга - но не 100%, всегда будут существовать зоны, свободные от точек. Да - это будет вида 0,00001% от общей площади - но будет.
Это как с интегралами. Площадь под кривой 3Х^2 от 1 до 2. Интеграл - легко видеть Х^3 и точное решение 8-1=7
А численными методами (суммируем площади прямоугольников под кривой) мы получим или 6,999999 или 7,000001 (смотря какие прямоугольники). Уменьшая толщину - мы повышаем точность = но точно 7 мы не получим.
PS
причем заданная точность вычислений (численные методы) будет работать в случае интеграла = нам нужна например точность 0,001%. уменьшаем толщину прямоугольников, пока разность между двумя вычислениями не составит требуемую величину.
а в случае случайного распределения точек мы не знаем исходной задачи, насколько критично наличие областей, в которых точек не будет вообще.
мы сначала переходим к пределу и получаем абсолютно точную (не на 0,0001%, а строго!) формулу для получения случайной равномерно-распределенной координаты точки внутри круга, а потом уже по этой точной формуле строим свои точки.
Да нет, там же математически показывается, что распределение у формулы точно такое же, как у предыдущей (с подстановкой в обратную CDF). Треугольник лишь нужен, чтобы додуматься до такой формулы.
Заполнить квадрат, обрезать лишнее до круга.
Так это тот же самый метод с отсечением (отбрасыванием), только вы шаг отсечения переносите в конец для всех точек. Плюс, сложнее понять, что мы закончили. Допустим, нам надо сгенерировать строго 100 точек внутри круга. Как мы поймём, что мы это сделали? Мы же не знаем, сколько из них в пределах круга, а сколько за пределами, пока не "обрежем".
Неплохо, для студента второго курса, изучающего тервер. Для продакшн имеет смысл использовать оба способа — вначале приближение целевой фигуры примитивом, для которого "легко" построить функцию равномерного распределения, затем проверка на попадание точки в фигуру.
а зачем нужен примитив?
проверка попадания x*x+y*y < R^2- в худшем случае- три умножения, сложение и сравнение.
примитив- определение сектора, в котором используем примитив (два сравнения- по xy), определение координаты внутри примитива, отображение на круг и снова проверка попадания с теми же накладными расходами. в чем профит? Если мне надо простой читаемый код, вызываемый сто раз в секунду- то используем равномерный квадрат с отбрасыванием- все просто и понятно, а если надо миллиарды раз в секунду такие точки строить- то я бы тогда больше не за отбрасывание переживал, а за синусы и косинусы- (или любые другие трансцендентные функции)- и снова вернулся бы к отбрасыванию- пересчитать рандом каждые 3 из 10 случаев быстрее и проще, чем два синуса в каждом случае.
Не в этой задаче, тут и так примитивом является или квадрат, или круг, и собственно вся статья посвящена способам сгенерировать равномерное распределение для примитива "круг". Я о более сложных задачах, где дана фигура неопределенной формы, но со списком условий на проверку принадлежности точки фигуре, и в ней нужно сделать равномерное распределение.
Забавно. Ровно эту же задачу я решал, делая для игры Arma 3 скрипт бомбардировки участка в определенном радиусе. Выбирал именно второй метод с полярными координатами. Помню, что заметил странную "равномерность" в распределения бомб, но списал это на "показалось".
Скорости выполнения получились примерно одинаковы, потому что весь пар уходит в вычисление тригонометрических функций. Стоит побороться за то, чтобы от них избавиться.
Навскидку вижу такой способ, который, правда, требует вычисления квадратного корня, что гораздо дешевле синуса, но тоже дороговато. Рассмотрим, как распределено значение координаты x. Ее плотность должна быть пропорциональнa длине хорды, которая образуется прямой, заданной уравнением x=const. То есть нужно вычислять отображение sqrt(1-x^2). Осталось найти распределение координаты y при заданном значении x: оно будет распределено равномерно по всей длине хорды. И нужен еще один вызов квадратного корня для нахождения длины хорды.
Хм, у меня попытка обратить интегральную функцию распределения F(x)=2/π∫√(1-x2)dx привела к трансцендентному уравнению y+sin y=C, так что не уверен, что эта идея сработает.
Да, действительно, так не получится, я забыл проинтегрировать плотность.
Вроде бы если сделать замену x = cos t, то из корня получится просто синус, делаем из него косинус по формулам приведения, заменяем обратно - и получается треугольная функция распределения. Так что способ действительно рабочий и без тригонометрии. Если я в уме нигде не ошибся :)
вместо того, чтобы преобразовать некое распределение, мы просто отклоняем все результаты, находящиеся вне нужного нам интервала.
Я не специалист, поэтому мне не очевидно: а характер распределения (в данном случае — равномерного) после отбрасывания части результатов — сохраняется?
Для равномерного распределения это как раз верно. Но неверно в общем случае, но там и распределения в среднем имеют бесконечные хвосты.
Будет чуть более очевидно, если представить, что мы сперва накидываем много точек на квадрат, а потом лишние убираем, и обратить внимание, что результат по сути тот же.
Если же хочется унылое доказательство, то можно считать вероятность попасть в элемент площади х. Это будет вероятность сразу попасть в него, плюс вероятность попасть в молоко и следом в х, плюс вероятность дважды в молоко и следом в х... Прогрессию складываем, выходит х нормированное на площадь круга.
Если хочется нормальное доказательство, то выйдет равномерно в силу симметрии: попадание в молоко не влияет на то, куда полетит следующая точка. Другими словами, в параграфе выше ответ зависит от площади х, но не зависит от того, где эта площадь.
Если помнить, что суммирование двух случайных величин эквивалентно свёртке их PDF, то код можно сделать ещё проще:
def sum_dist():
theta = random() * pi
r = random() - random()
return r * cos(theta), r * sin(theta)— так как вычисление random() - random () даёт случайную величину с треугольным распределением, симметричным относительно 0.
Мне на ум сразу пришло одно наиболее востребованное применение - онлайн казино, рулетка)
Есть еще такой вариант, не проверял на скорость но по моему он может быть даже быстрее чем Rejection sampling. Не знаю как назвать этот вариант, она просто изменяет область значений y основываясь на значении x.
def nosin_dist(n):
x = random() * 2 - 1;
y = (random() * 2 - 1)*sqrt(1 - x^2);
return x, y;Такой способ не даст равномерное распределение, плотность в окрестности точек (-1,0) и (1,0) будет заметно выше.
Правильно, я пытался найти функцию плотности вероятности, только моих мат. знаний оказалось недостаточно.
Вот тут уже пытались. Это совсем не тривиально.
Тут та же проблема что и в наивном методе в полярных координатах — не равномерное распределение.
в форме твоего видео для 3Blue1Brown было интересней.
Значительно быстрее остальных оказалась выборка с отклонением (rejection sampling)
Не совсем так. В среднем, да.
А если считать для систем реального времени, то мы должны закладывать самую долгий вариант расчёта. А он тут бесконечность.
Лучший способ выбора случайной точки в круге