Engee 文档
Notebook

数学统计

在本示例中,我们将演示如何解决数理统计领域的简单问题--获取一系列随机数、建立离散分布模型、构建正态分布随机变量的概率密度以及卡尔曼滤波示例。

首先,执行准备单元格以设置图形的绘制。

In [ ]:
Pkg.add(["Distributions", "LinearAlgebra", "StatsBase", "Measures"])
In [ ]:
using LaTeXStrings;
using Plots

gr()
#plotlyjs(); # Для интерактивных графиков
Out[0]:
Plots.GRBackend()

实现随机变量

任务: 从方差为 1、数学期望等于-2 的随机变量的正态分布中抽取 100 个点。

In [ ]:
n = 100;
ϵ = randn(n) .- 2;
plot(1:n, ϵ, label=L"\mathcal{N}(-2, 1)")
Out[0]:

离散分布:骰子

任务: 比较一个随机变量的两种分布,该变量描述了投掷一个六角骰子的结果,是通过蒙特卡罗方法和库Distributions 得到的。

In [ ]:
using Distributions, StatsBase, Plots;

faces, N = 1:6, 10^4
mcEstimate = counts(rand(faces,N), faces)/N

a = Distributions.DiscreteUniform(1,6)

plot(faces, mcEstimate, line=:stem, marker=:circle, c=:blue,
    ms=10, msw=0, lw=8, label="Монте-Карло", xticks=1:12, )

plot!([i for i in faces], [pdf(a, x) for x in faces], line=:stem, label="Плотность",
    lw=5, lc=:skyblue,
    xlabel="Выпавшее число", ylabel="Вероятность", ylims=(0,0.22))
Out[0]:

正态随机变量的概率密度

任务: 绘制数学期望为 1.5、方差等于 3.1 的分布的概率密度函数。

In [ ]:
using Distributions, Plots, LaTeXStrings
d = Normal(1.5, 3.1)

x = sort!(rand(d,1000))
y = pdf.(d, x)

plot(x, y, label=L"\frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right)")
ylabel!(L"f(x; \mu, \sigma)")
Out[0]:

卡尔曼滤波的一维示例 [1,p.343]

许多物理系统都可以用下面的模型来描述:

$$\begin{matrix} X(t+1) &= &AX(t) + \xi(t)\\ Y(t) &= &CX(t) + \zeta(t)\end{matrix}$$

卡尔曼滤波器中的增益矩阵可以通过迭代计算来改善滤波结果。在这个例子中,所有变量都是标量。我们考虑$a \in (0,1)$ ,模型会随着$a$ 在每个时间瞬间的强度趋于 0。

任务: 演示卡尔曼滤波跟踪单变量自回归过程的能力。

In [ ]:
using Distributions, LinearAlgebra, Random, Measures, Plots;
Random.seed!(1)

a, varXi, varZeta, = 0.8, 0.36, 1.0
X0, spikeTime, spikeSize = 10.0, 50, -20.0
Tsmall, Tlarge, warmTime = 100, 10^6, 10^4
kKalman = 0.3

function luenbergerTrack(k, T, spikeTime = Inf)
    X, Xhat = X0, 0.0
    xTraj, xHatTraj, yTraj = [X], [Xhat,Xhat], [X0]

    for t in 1:T-1
        X = a*X + rand(Normal(0,sqrt(varXi)))
        Y = X + rand(Normal(0,sqrt(varZeta)))
        Xhat =a*Xhat - k*(Xhat - Y)
        
        push!(xTraj,X)
        push!(xHatTraj,Xhat)
        push!(yTraj,Y)
    
        if t == spikeTime
            X += spikeSize
        end
    end
    deleteat!(xHatTraj,length(xHatTraj))
    xTraj, xHatTraj, yTraj
end

smallTraj, smallHat, smallY = luenbergerTrack(kKalman, Tsmall, spikeTime)

p1 = scatter(smallTraj, c=:blue, ms=3, msw=0, label="Траектория системы")
p1 = scatter!(smallY, c = :black, ms=3, msw=0, label="Измерения")
p1 = scatter!(smallHat,c=:red, ms=3, msw=0, label="Kalman filter tracking",
xlabel = "Время", ylabel = "Температура", xlims=(0, Tsmall))

kRange = 0.2:0.005:0.4
errs = []
for k in kRange
    xTraj, xHatTraj, _ = luenbergerTrack(k, Tlarge)
    mse = norm(xTraj[warmTime:end] - xHatTraj[warmTime:end])^2/(Tlarge-warmTime)
    push!(errs, mse)
end

analyticErr(k) = (varXi + k^2*varZeta) / (1-(a-k)^2)

p2 = scatter(kRange,errs, c=:black, ms=3, msw=0,
xlabel="k", ylabel="СКО", label = "Monte Carlo")
p2 = plot!(kRange,analyticErr.(kRange), c = :red,
xlabel="k", ylabel="СКО", label = "Аналитик", ylim =(0.58,0.64))

plot(p1, p2, size=(1000,400), margin = 5mm)
Out[0]: