Документация Engee
Notebook

Определение параметров распределения по частичной выборке

Открыть пример в Engee

В этом примере мы разделим выборку на две части и аппроксимируем каждую половину отдельной функцией плотности вероятности, воспользовавшись методом наименьших квадратов для нахождения параметров каждой функции.

Создание выборки

Нашей выборкой будет модель броуновского процесса на отрезке [-20,20]. Переход частицы происходит с равной вероятностью влево и вправо, длина пробега определяется по экспоненциальному распределению с заданным единичным средним.

Получим набор реализаций этого процесса, запустив 100000 траекторий из центра области.

In [ ]: 
a = -20
b = 20
λ = 1

using Random
Random.seed!(2)

Ntraj = 100_000  # Количество рассчитываемых траекторий
Nc = []          # Набор длин траекторий
NtrajList = []   # Список координат траектории частиц

for i in 1:Ntraj
    x = 0
    k = 0
    if length(NtrajList) < 200 push!( NtrajList, [] ); end
    while a < x < b
        r1,r2 = rand(2)
        dx = -log(r2) / λ
        if( r1 < 0.5 ) x -= dx; else x += dx; end;
        k += 1;
        # Сохраним первые 200 траекторий в память
        if length(NtrajList) < 200 push!( NtrajList[end], x ); end
    end
    Nc = push!( Nc, k )
end

Nc = Float64.( Nc );

Нам пришлось перевести вектор Nc в тип Float64, без чего мы часто получали бы предупреждения о том, что в векторе могут находиться данные типа Nothing или Missing, даже если фактически их там нет.

Выведем первые пару сотен траекторий.

In [ ]: 
gr()

plot()
for traj in NtrajList
    plot!( 1:length(traj), traj, c=1 )
end

using Measures
plot!( legend=:false, size=(1000,400), xlabel="Шаги симуляции", ylabel="Положение частицы", left_margin=10mm, bottom_margin=10mm )
Out[0]:

Мы видим, что некоторые частицы очень долго находились в контролируемом пространстве (до 800 тактов и более), но в основном они покинули его в течение первой сотни шагов от запуска симуляции.

Найдем статистические параметры этой выборки:

In [ ]: 
using Statistics
Nave = mean( Nc ); # Среднее значение
Nd = std( Nc );    # Дисперсия

histogram( Nc, xlims=(0,1000), normalize=true, label="Данные" )
vline!( [Nave], lw=4, lc=2, label = "Среднее" )
vline!( [Nave+Nd], lw=3, lc=3, label=:none )
vline!( [Nave-Nd], lw=3, lc=3, label=:none )
Out[0]:

Чтобы исследовать эту выборку, среднего арифметического значения недостаточно. Нам нужно угадать, какой функции подчинено распределение результатов наблюдаемого процесса. Для начала предположим, что перед нами логнормальное распределение.

In [ ]: 
using Distributions

dist0 = fit( LogNormal, Nc )
x1 = 10:1000

histogram( Nc, xlims=(0,1000), normalize=true, label="Данные" )
plot!( x1, pdf.( dist0, x1 ), lw=4, lc=2, label="Логнормальное распределение" )
Out[0]:

Правый хвост гистограммы по своему виду заметно отличается от кривой найденного логнормального распределения.

Разделим данные на две выборки

Попробуем изучить две выборки отдельно одну от другой. Для этого разделим выборку на две половины (до и после максимума) и найдем параметры двух функций распределения. Прежде всего, где находится максимум распределения?

In [ ]: 
xx = range( 10, 1000, step=8 )        # Интервалы гистограммы
xc = xx[1:end-1] .+ (xx[2]-xx[1])/2   # Центры интервалов

#dist_h = fit( Histogram, Nc, xx ).weights; # Максимум гистограммы дает не лучшую оценку из-за выбросов,
dist_h = pdf.( dist0, x1 )                  # поэтому возьмем максимум логнормальной функции

m,i = findmax( dist_h )
Nc_mid = x1[i]
Out[0]:
88

Отделим две подвыборки:

In [ ]: 
Nc1 = Nc[ Nc .<= Nc_mid ];
Nc2 = Nc[ Nc .> Nc_mid ];

Аппроксимируем эти две выборки независимо одну то другой. Построим по ним гистограммы. Обратите внимание на аргумент bins – если его не указывать, система выберет разный размер интервалов для каждой выборки, а поскольку выборки у нас разного размера, график плотности будет ненормированным относительно размера интервала.

In [ ]: 
histogram( Nc1 )
histogram!( Nc2, size=(600,150) )
Out[0]:
In [ ]: 
histogram( Nc1, bins=1:10:Nc_mid, label="Левая половина выборки" )
histogram!( Nc2, bins=Nc_mid+1:10:1000, label="Правая половина выборки" )
Out[0]:

Теперь у нас есть веса каждого интервала, полученные после построения гистограмм.

In [ ]: 
using StatsBase

# Границы интервалов для гистограмм
x1 = 10:Nc_mid
x2 = Nc_mid+1:1000

# Центры интервалов
c1 = x1[1:end-1] .+ (x1[2]-x1[1])/2
c2 = x2[1:end-1] .+ (x2[2]-x2[1])/2

# Значения гистограммы (третий параметр – границы интервалов)
w1 = fit( Histogram, Nc1, x1 ).weights;
w2 = fit( Histogram, Nc2, x2 ).weights;

# Нормализованные гистограммы (важно нормализовать используя всю длину Nc, а не длину отдельных хвостов)
w1n = w1 ./ length(Nc);
w2n = w2 ./ length(Nc);

plot( c1, w1n, label="Левая половина выборки" )
plot!( c2, w2n, label="Правая половина выборки" )
Out[0]:

Но каждый из них представляет лишь часть данных для построения гистограммы. По этим данным нельзя полноценно найти параметры распределения, поскольку перед нами ограниченная выборка.

In [ ]: 
plot( c1, w1n, label="Левая половина выборки" )
plot!( c2, w2n, label="Правая половина выборки", size=(600,180) )

dist1 = fit( LogNormal, Nc1 )
dist2 = fit( Exponential, Nc2 )

plot!( c1, pdf.( dist1, c1 ), lw=4, label="Логнормальное распределение" )
plot!( c2, pdf.( dist2, c2 ), lw=4, label="Экспоненциальное распределение" )
Out[0]:

Но зато мы можем найти параметры распределения, если задать его как функцию, которая должна пройти через представленные точки.

Параметры левой выборки

Представим себе, что левая часть распределения следует логнормальному закону. Это распределение имеет следующую формулу плотности вероятности:

$$P(x) = \exp \left ( - \left [ \frac{ \ln (x) - \mu } {\sigma}\right ]^2 \bigg/ 2 \right ) \Bigg / (x \sigma \sqrt{2 \pi})$$

Эта формула имеет два параметра – $\mu$ и $\sigma$, которые у нас будут образовывать элементы вектора p (p[1] и p[2]).

Найдем параметры этой кривой с условием, что она должна максимально близко проходить через имеющиеся точки.

In [ ]: 
# Впишем функцию средствами подбора параметров функций (а не распределений)
using LsqFit

# Эти функции должны быть написаны для векторных аргументов, поэтому @. вначале (чтобы все операции были .*, .+ и т.д.)
@. lognormDist( x, p ) = (1 / (x * p[2] * sqrt(2π) )) * exp(-(log(x) - p[1])^2 / (2p[2]^2));

p0 = [1.0, 1.0];
af1 = curve_fit( lognormDist, c1, w1n, p0 );

af1.param
Out[0]:
2-element Vector{Float64}:
 5.0813235614097465
 0.8009541756702706

Мы подобрали параметры распределения, которые дают нам плотность вероятности, очень похожую на ту, которую нам дает встроенная функция LogNormal.

In [ ]: 
plot( c1, lognormDist(c1, af1.param), label="Наша функция LogNorm", legend=:topleft )
plot!( c1, pdf.(LogNormal(af1.param[1], af1.param[2]), c1), label="Штатная функция для плотности Lognorm" )
plot!( c1, w1n, label="Гистограмма по которой мы строили распределение", legend=:bottomright )
Out[0]:

Параметры правой выборки

Попробуем вписать точки "правой" половины выборки в экспоненциальную функцию плотности вероятности:

$$P(x) = \lambda_1 e^{-\lambda_2 x}$$

Иногда это распределение строится через однопараметрическую функцию, но в нашем случае результаты получаются лучше, если $\lambda_1 \neq \lambda_2$.

In [ ]: 
@. expDist( x, p ) = p[1] * exp( -p[2]*x );

p0 = [ 0.1, 0.1]; # Все очень сильно зависит от начального параметра. При значениях выше 0.5 модель не сходится
af2 = curve_fit( expDist, c2, w2n, p0 );

af2.param
Out[0]:
2-element Vector{Float64}:
 0.00685261267119607
 0.005432207598447344

Поскольку мы пользуемся итерационным методом для аппроксимации выборки при помощи заданной функции, нужно указывать исходную аппроксимацию решения p_0, от выбора которой часто зависит, сойдется или нет наш алгоритм.

In [ ]: 
plot( c2, expDist(c2, af2.param), label="Экспоненциальное распределение" )
plot!( c2, w2n, label="Гистограмма, по которой мы искали параметры" )
Out[0]:

Совместим обе части

Если нанести обе функции на один график, мы получим довольно хорошее соответствие с исходной гистограммой, в чем можем убедиться при помощи пакета HypothesisTests.

In [ ]: 
stephist( Nc, xlims=(10,1000), bins = 10:5:1000, normalize=true, label="Исходные данные", fill=true, fillalpha=0.3, color=:black, fillcolor=1 )
plot!( [c1; c2], pdf.( dist0, [c1; c2] ), lw=3, lc=2, ls=:dash, label="Логнормальное распределение" )
plot!( [c1; c2], [lognormDist(c1, af1.param); expDist(c2, af2.param)], lw=3, lc=:black, ls=:dash, label="Найденная плотность распределения" )
Out[0]:

Оценим качество этого решения

Мы фактически вписали кусочно-линейную функцию в выборку. Теперь оценим остаточную ошибку – как визуально, так и при помощи критерия Стьюдента.

In [ ]: 
using HypothesisTests

yy1 = pdf.( dist0, [c1; c2] )
yy2 = [ lognormDist(c1, af1.param); expDist(c2, af2.param) ]
yy3 = [w1n; w2n]

plot(
    plot( yy1 .- yy3, size=(600,250), lc=2,
        label="Lognorm (pvalue = $(round(pvalue(OneSampleTTest(yy1 .- yy3)),digits=3)))" ),
    plot( yy2 .- yy3, size=(600,250), lc=1,
        label="Lognorm + Exponential (pvalue=$(round(pvalue(OneSampleTTest(yy2 .- yy3)),digits=3)))" ),
    layout=(2,1)
)
Out[0]:

Второй график ошибки выглядит более однородно.

$p$-значение показывает, что мы учли больше систематических эффектов, и что ошибка между выборкой и составленной нами функцией ближе к нормально распределенной случайной величине, чем ошибка между выборкой и функцией распределения LogNorm.

Заключение

Мы решили прикладную задачу, для которой нам потребовалось вписать несколько функций в выборки. В результате мы построили сложную функцию распределения для описания результатов физического эксперимента.