Распределение случайной величины в разных системах координат¶
В данном примере будут рассмотрены следующие темы:
- генерация случайных чисел и массивов в Julia
- использование распределений вероятности
- важность правильного выбора распределения случайной величины в зависимости от постановки задачи.
Генерация случайных чисел по умолчанию¶
функция rand()¶
Для того чтобы сгенерировать случайное число, можно без подключения библиотек воспользоваться функцией rand() . Она вернёт случайное вещественное число в диапазоне [0,1)
rand()
Если нужно создать матрицу $A^{i\times j \times k}$ определённого размера с элементами $a_{ijk}\in[0,1)$, то воспользуемся rand(i,j,k)
rand(2,4,3)
Если же в rand() передадать кортеж:
rand((2,4,3)),
то (2,4,3) будет множеством элементов, из которых выбирается случайный элемент
rand((2,4,3))
rand(('a',3,"string",1+1im),2,3)
Также в качестве множества элементов можно передавать словарь:
rand(Dict("pi" => 3.14, "e" => 2.71), 3)
Визуальный анализ распределения¶
Для того чтобы визуально оценить, по какому закону генерируются числа, воспользуемся библиотекой Plots. Для этого подключим её с помощью using Plots и после этого воспользуемся функцией histogram. Прописав gr() - отключим возможность интерактивного взаимодействия с графиками.
using Plots # подключаем библиотеку
gr() # отключаем интерактивное взаимодействие с графиком, т.к. точек много
histogram(rand(50000)) # строим гистограму
title!("rand(50000)") # подпись
Как видно из гистограмы, распределение случайных чисел, сгенерированных при помощи rand(), является равномерным. Т.е. функция rand() генерирует величины со стандартным непрерывным равномерным распределением (СНРР)
Нормальное распределение¶
Также, без подключения других библиотек, в Julia можно генерировать случайные числа с нормальным распределением.
Функция randn() генерирует нормальное распределение чисел с математическим ожиданием 0 и среднеквадратичным отклонением 1. (т.е. диапазон значений будет $(-\infty,\infty$)
histogram(randn(50000),bims=0.005) # bims - ширина стобца, в который попадает точка
title!("randn(50000)")
Важность выбора правильного распределения случайной величины¶
Рассмотрим такую задачу:
Создать набор точек внутри круга радиуса R=1 и с центром в (0,0), декартовы координаты которых равномерно распределены по осям абсцисс и ординат. Для этого
- Создадим набор точек, попадающих в квадрат с диагональю , проходящей через вершины
(-1,1), (1,-1) - Выберем те точки, которые которые попадают в единичный круг.
Шаг 1.
Для решения этого шага, необходимо помнить, что мы хотим получить равномерно распределённые значения X в диапазоне [-1,1]. А rand() - СНРР (т.е. [0,1)).
СНРР обладает свойством масштабирования: Если случайная величина $X \sim U[0,1]$ и $Y = a+(b-a)X$ , то $Y \sim U[\min(a,b),\max(a,b)]$
В нашем случае $a=-1$ и $b=1$ $\Rightarrow$
POINTS_NUMBER = 3000
big_square_points = -1 .+ 2. .* rand(POINTS_NUMBER, 2) # воспользовались свойством СНРР
Для того чтобы при построении графиков каждый раз не использовать одни и те же ключевые слова можно добавить PlotThemes и прописать необходимые слова "по умолчанию".
Pkg.add("PlotThemes")
theme(:default, aspect_ratio=1, label = false, markersize=1)
scatter(big_square_points[:,1],big_square_points[:,2])
Шаг 2
Будем проходить по всем точками нашего квадрата и добавлять в массив те координаты, которые попадают в единичный круг.
circle_points_x = Float64[]
circle_points_y = Float64[]
R = 1
for (x, y) in eachrow(big_square_points) # (в каждой строке big_square_points 2 элемента: x и y)
if x^2 + y^2 < R^2
push!(circle_points_x, x)
push!(circle_points_y, y)
end
end
scatter(circle_points_x,circle_points_y)
Для решение нашей задачи более естественным подходом было бы использовать не декартовы координаты, а полярные.
Использование полярных координат¶
Тогда эта задача сведётся к одному, а не двум шагам:
Шаг 1 Выбрать $\rho \in [0,1]$ и $\alpha \in [0, 2\pi]$
ρ = rand(POINTS_NUMBER)
α = 2π * rand(POINTS_NUMBER)
polar_original = scatter(ρ .* cos.(α), ρ .* sin.(α), aspect_ratio=1,markersize=1) # перешли обратно в декартовы
Также можно сразу посмотреть наши результаты, без перехода в декартовы координаты. Для этого предварительно отключим aspect_ratio=1 и укажем proj=:polar - полярное отображение.
theme(:default)
scatter(α, ρ , proj=:polar, markersize=1)
Можем заметить, что плотность точек возле начала координат выше, чем возле окружности.
Для наглядности вырежем из окружности вписанный квадрат.
Для этого будем использовать функцию zip():
for i in zip(1:5,'a':'e',["🟥","🟩","🟦"])
println(i)
end
# 4,5 и 'd','e' не напечатаются, так как коллекция цветов закончилась на 3 шаге
Визуализация плотности точек¶
small_square_points_x = Float64[]
small_square_points_y = Float64[]
for (x, y) in zip(circle_points_x, circle_points_y)
if abs(x) ≤ √2 / 2 && abs(y) ≤ √2 / 2
push!(small_square_points_x, x)
push!(small_square_points_y, y)
end
end
small_square_points_ρ = Float64[]
small_square_points_α = Float64[]
for (x, y) in zip(ρ .* cos.(α), ρ .* sin.(α))
if abs(x) ≤ √2 / 2 && abs(y) ≤ √2 / 2
push!(small_square_points_ρ, x)
push!(small_square_points_α, y)
end
end
Чтобы показать сразу несколько графиков на одном полотне необходимо использовать layout и указать структуру их расположения (в нашем случае сетка $2 \times 2$. При помощи histogram2d сможем увидеть, в какой области большая плотность точек.
p1 = scatter( small_square_points_x, small_square_points_y,
markersize=1, legend = false, aspect_ratio=1)
p2 = scatter( small_square_points_ρ, small_square_points_α,
markersize=1, legend = false, aspect_ratio=1)
p3 = histogram2d(small_square_points_x,small_square_points_y,
title="плотность в декартовых координатах", aspect_ratio=1)
p4 = histogram2d(small_square_points_ρ,small_square_points_α,
title="плотность в полярных координатах",aspect_ratio=1)
plot(p1,p2,p3,p4,layout=(2,2),size=(950,950))
Функция repeat()¶
Рассмотрим функцию repeat()
repeat([1:3;],2)
repeat([1:3;],1,2)
repeat([1:3;],2,3)
repeat([1:3;],inner=2)
Сравнение полярных и декартовых координат¶
Чтобы понять, почему применение полярных координат и случайное равномерное распределение по $\rho$ и $\alpha$ сработали таким, а не, возможно, ожидаемым образом, рассмотрим следующую картину. Возьмём не случайное, а фиксированное равномерное распределение:
lngt = 12
a = repeat(range(0, stop=2pi,length=lngt),1,lngt)
r = repeat(range(0, stop=1, length=lngt),1,lngt)'
scatter(r .* cos.(a), r .* sin.(a),aspect_ratio=1, markersize=3)
x = repeat(range(0, stop=1, length=lngt),1,lngt)'
y = repeat(range(0, stop=1,length=lngt),lngt)
scatter(x, y,aspect_ratio=1)
Как видим, плотность распределения точек в декартовых координатах является постоянной, а в полярных координатах плотность ваше возле начала координат.
Чтобы решить задачу нам нужно использовать такое распределение, которое будет с большей вероятностью ставить точку ближе к концу радиус-вектора, чем к началу. Предположим, что такое распределение должно быть линейным (чем ближе к концу, тем вероятнее). Оно называется треугольным.
Треугольное распределение¶
Чтобы воспользоваться треугольным распределением, подключим библиотеку Distributions. Данная библиотека имеет большой набор различных распределний (в частности, одномерных). Ознакомиться с этим набором можно в документации.
using Pkg
Pkg.add("Distributions")
using Distributions
Теперь создадим набор точек с треугольным распределением. Для этого в качестве первого аргумента функции rand() передадим TriangularDist(a,b,c) с параметрами a=0,b=1,c=1.
histogram(rand(TriangularDist(0,1,1),POINTS_NUMBER))
Убедившись, что треугольное распределение работает так, как ожидалось,проверим, что его применение решит нашу задачу:
ρ_triang = rand(TriangularDist(0,1,1), POINTS_NUMBER)
α_triang = rand(Uniform(0,2π), POINTS_NUMBER) # Uniform(0,2π) - равномерное распределение от 0 до 2π
polar_triang = scatter(ρ_triang .* cos.(α_triang), ρ_triang .* sin.(α_triang),aspect_ratio=1, markersize=1)
plot(polar_original, polar_triang,layout=(1,2))
В качестве интересного факта, покажем, как без подключения библиотеки Distributions можно было решить эту задачу, добавив всего два символа .√ :
ρ_sqrt = .√rand(POINTS_NUMBER)
α_sqrt = 2π * rand(POINTS_NUMBER)
scatter(α_sqrt, ρ_sqrt, proj=:polar, markersize=1)
Это справедливо из-за того, что:
Используя генератор равномерного распределения, можно получить генератор треугольного распределения: $$X = \begin{cases} a + \sqrt{U(b-a)(c-a)} & \text{ для } 0 < U < F(c) \\ & \\ b - \sqrt{(1-U)(b-a)(b-c)} & \text{ для } F(c) \le U < 1, \end{cases}$$ где $U$ - равномерно распределённое значение, а $F(c) = (c-a)/(b-a)$.
Подставив a=0,b=1,c=1, как раз и получим .√rand()
Заключение¶
В данном примере было показано, что язык Julia имеет как удобные встроенные функции, так и специализированные пакеты, позволяющие использовать специфические для области функции, взаимодействуя со стандартными.
Также были рассмотрены способы генерации массивов (repeat) и итерации по ним (zip и eachrow)