Математическая статистика
В этом примере мы продемонстриурем, как решать простые задачи из области математической статистики – получение ряда случайных чисел, моделирование дискретного распределения, построение плотности вероятности нормально распределенной случайной величины и пример фильтрации Калмана.
Сперва выполните подготовительные ячейки для настройки отрисовки графиков.
using LaTeXStrings;
using Plots
gr()
#plotlyjs(); # Для интерактивных графиковPlots.GRBackend()Реализация случайной величины
Задача: Отобразите сто точек из нормального распределения случайной величины с дисперсией 1 и математическим ожиданием, равным -2.
n = 100;
ϵ = randn(n) .- 2;
plot(1:n, ϵ, label=L"\mathcal{N}(-2, 1)")
Дискретное распределение: игральный кубик
Задача: Сравнить два распределения случайной величины, которая описывает результат броска одного шестигранного кубика, полученные методом Монте-Карло и через библиотеку Distributions.
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))
Плотность вероятности нормальной случайной величины
Задача: Построить график плотности вероятности распределения с математическим ожиданием 1.5 и дисперсией, равной 3.1.
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)")
Одномерный пример фильтрации Калмана [1,c.343]
Многие системы физической природы можно описать при помощи модели следующей формы:
Матрицы усиления в фильтре Калмана можно рассчитывать итерационно, чтобы улучшать результаты фильтрации. В этом примере все переменные являются скалярными величинами. Мы считаем, что и модель стремится к 0 с силой в каждый момент времени.
Задача: Продемонстрируйте способность фильтра Калмана отслеживать одномерный авторегрессионный процесс.
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)