Engee 文档
Notebook

随机变量在不同坐标系中的分布

本例将涉及以下主题：

  • 在 Julia 中生成随机数和数组
  • 使用概率分布
  • 根据问题陈述选择正确的随机变量分布的重要性。

默认随机数生成

函数rand()

为了生成随机数，您可以使用rand() 函数，而无需连接库。它将返回一个范围在[0,1)

In [ ]: 
Pkg.add(["Distributions"])
In [ ]: 
rand()
Out[0]:
0.8968475595278809

如果我们需要创建一个具有一定大小的矩阵$A^{i\times j \times k}$ ，其元素为$a_{ijk}\in[0,1)$ ，我们将使用rand(i,j,k)

In [ ]: 
rand(2,4,3)
Out[0]:
2×4×3 Array{Float64, 3}:
[:, :, 1] =
 0.105281  0.553249  0.571817  0.068993
 0.463789  0.510086  0.468581  0.599769

[:, :, 2] =
 0.845341  0.959642  0.488022  0.959579
 0.409156  0.49006   0.548231  0.899735

[:, :, 3] =
 0.957553  0.136052  0.885431  0.284849
 0.826566  0.92279   0.659209  0.384685

如果我们向rand() 传递一个元组：

rand((2,4,3)),

那么(2,4,3) 将是一个元素集合，从中随机抽取一个元素

In [ ]: 
rand((2,4,3))
Out[0]:
4
In [ ]: 
rand(('a',3,"string",1+1im),2,3)
Out[0]:
2×3 Matrix{Any}:
   "string"  'a'    "string"
 1+1im       'a'  1+1im

您也可以传递一个字典作为元素集：

In [ ]: 
rand(Dict("pi" => 3.14, "e" => 2.71), 3)
Out[0]:
3-element Vector{Pair{String, Float64}}:
 "e" => 2.71
 "e" => 2.71
 "e" => 2.71

可视化分布分析

为了直观地评估数字是根据哪种规律生成的，让我们使用Plots 库。为此，我们使用using Plots 将其连接起来，然后使用函数histogram 。通过指定gr() ，我们将禁用与图表交互的可能性。

In [ ]: 
using Plots             # подключаем библиотеку
gr()                    # отключаем интерактивное взаимодействие с графиком, т.к. точек много
histogram(rand(50000))  # строим гистограму 
title!("rand(50000)")   # подпись
Out[0]:

从直方图中我们可以看到，使用 rand() 生成的随机数的分布是均匀的。也就是说，rand() 函数生成的数值是标准连续均匀分布（SNRD）。

正态分布

在 Julia 中也可以生成正态分布随机数，而无需连接其他库。

randn()函数生成正态分布的数字，其数学期望为0 ，标准偏差为1 。(即数值范围为$(-\infty,\infty$)

In [ ]: 
histogram(randn(50000),bims=0.005) # bims - ширина стобца, в который попадает точка
title!("randn(50000)")
Out[0]:

选择正确的随机变量分布的重要性

让我们考虑这样一个问题：

在一个半径为R=1 、圆心位于(0,0) 的圆内创建一组点，这些点的笛卡尔坐标沿 abscissa 轴和 ordinate 轴均匀分布。为此 1.在对角线 , 经过顶点的正方形中创建一组点。(-1,1), (1,-1) 2. 选择落在单位圆内的点。

步骤 1.

要解决这一步问题，请记住我们希望X 的值在[-1,1] 范围内均匀分布。而rand() 是 SNRR（即[0,1) ）。

SNRR 具有 缩放特性： 如果随机变量$X \sim U[0,1]$ 和$Y = a+(b-a)X$ ，则$Y \sim U[\min(a,b),\max(a,b)]$

在我们的情况下，$a=-1$ 和$b=1$ $\Rightarrow$

In [ ]: 
POINTS_NUMBER = 3000
big_square_points = -1 .+ 2. .* rand(POINTS_NUMBER, 2)  # воспользовались свойством СНРР
Out[0]:
3000×2 Matrix{Float64}:
  0.0529451    -0.378398
 -0.926168     -0.595867
 -0.506066      0.586659
 -0.349772     -0.501738
 -0.199344     -0.742828
  0.869969      0.759997
  0.0879778    -0.111934
 -0.601932     -0.955655
  0.678965      0.0709128
 -0.873571     -0.671576
  0.161069      0.667479
  0.32568       0.848502
  0.14866      -0.0551039
  ⋮            
 -0.403644     -0.839549
 -0.103682     -0.11301
 -0.958197     -0.63964
  0.0131706     0.507728
  0.0394808     0.353285
 -0.948482     -0.794558
 -0.870387     -0.195976
  0.980847      0.716624
 -0.987589      0.0785435
  0.897117     -0.967298
  0.000805525  -0.625087
 -0.112197     -0.540335

为了避免在构建图表时每次都使用相同的关键字，可以添加PlotThemes 并 "默认 "指定必要的字词。

In [ ]: 
Pkg.add("PlotThemes")
theme(:default, aspect_ratio=1, label = false, markersize=1)
scatter(big_square_points[:,1],big_square_points[:,2])  

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
Out[0]:

第 2 步

我们将遍历正方形的所有点，并将落在单位圆内的坐标添加到数组中。

In [ ]: 
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)
Out[0]:

要解决我们的问题，更自然的方法是使用极坐标而不是笛卡尔坐标。

使用极坐标

那么这个问题就可以简化为一个步骤，而不是两个步骤：

**第 1 步 选择$\rho \in [0,1]$ 和$\alpha \in [0, 2\pi]$

In [ ]: 
ρ = rand(POINTS_NUMBER)
α = 2π * rand(POINTS_NUMBER)
polar_original = scatter(ρ .* cos.(α), ρ .* sin.(α), aspect_ratio=1,markersize=1)  # перешли обратно в декартовы
Out[0]:

也可以不切换到直角坐标，直接查看结果。为此，我们禁用aspect_ratio=1 并指定proj=:polar - 极坐标显示。

In [ ]: 
theme(:default)
scatter(α, ρ , proj=:polar, markersize=1)
Out[0]:

我们可以发现，原点附近的点密度高于圆附近的点密度。

为了说明问题，让我们从圆中切出一个内切正方形。

为此，我们将使用函数zip() ：

In [ ]: 
for i in zip(1:5,'a':'e',["🟥","🟩","🟦"])
    println(i)
end
# 4,5 и 'd','e' не напечатаются, так как коллекция цветов закончилась на 3 шаге

(1, 'a', "🟥")
(2, 'b', "🟩")
(3, 'c', "🟦")

点密度可视化

In [ ]: 
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
In [ ]: 
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 的帮助下，我们可以看到哪个区域的点密度较高。

In [ ]: 
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))
Out[0]:

函数repeat()

考虑函数repeat()

In [ ]: 
repeat([1:3;],2)
Out[0]:
6-element Vector{Int64}:
 1
 2
 3
 1
 2
 3
In [ ]: 
repeat([1:3;],1,2)
Out[0]:
3×2 Matrix{Int64}:
 1  1
 2  2
 3  3
In [ ]: 
repeat([1:3;],2,3)
Out[0]:
6×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3
 1  1  1
 2  2  2
 3  3  3
In [ ]: 
repeat([1:3;],inner=2)
Out[0]:
6-element Vector{Int64}:
 1
 1
 2
 2
 3
 3

极坐标与直角坐标的比较

要理解为什么在$\rho$ 和$\alpha$ 上应用极坐标和随机均匀分布会有这种效果，而不是预期的效果，请看下图。 让我们采用固定的均匀分布，而不是随机分布：

In [ ]: 
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)
Out[0]:
In [ ]: 
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)
Out[0]:

我们可以看到，在直角坐标中，点的分布密度是恒定的，而在极坐标中，你的密度在原点附近。

为了解决这个问题，我们需要使用这样一种分布，它更有可能使一个点靠近半径向量的末端而不是起点。假设这种分布应该是线性的（越靠近末端，可能性越大）。这就是三角形分布。

三角形分布

为了利用 三角形 分布，让我们插入Distributions 库。这个库有大量不同的分布（尤其是一维分布）。这组分布可以在 documentation 中找到。

In [ ]: 
using Pkg   
Pkg.add("Distributions")
using Distributions

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

现在我们来创建一组三角形分布的点。为此，作为函数rand() 的第一个参数，传递TriangularDist(a,b,c) 及参数a=0,b=1,c=1 。 triangular_distribution_pmf.png

In [ ]: 
histogram(rand(TriangularDist(0,1,1),POINTS_NUMBER))
Out[0]:

在验证了三角形分布符合预期之后，让我们来检验一下它的应用是否能解决我们的问题：

In [ ]: 
ρ_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))
Out[0]:

一个有趣的事实是，我们展示了如何在不连接库Distributions 的情况下，只添加两个符号.√ 就能解决这个问题：

In [ ]: 
ρ_sqrt = .√rand(POINTS_NUMBER)
α_sqrt = 2π * rand(POINTS_NUMBER)
scatter(α_sqrt, ρ_sqrt, proj=:polar, markersize=1)
Out[0]:

这是因为

利用均匀分布发生器，我们可以得到三角分布发生器： $$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 ）。