Engee 文档
Notebook

确定部分样本的分布参数

在这个例子中，我们将把样本分成两部分，并用一个单独的概率密度函数近似每一半，使用最小二乘法找到每个函数的参数。

创建选择

我们的样本将是区间[-20,20]中布朗运动过程的模型。 粒子转变以相等的概率向左和向右发生，路径长度由具有给定单位平均值的指数分布确定。

我们将通过从区域中心发射100,000个轨迹来获得这一过程的一组实现。

In [ ]: 
Pkg.add(["Distributions", "Statistics", "Random",
         "StatsBase", "Measures", "HypothesisTests", "LsqFit"])
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 (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]:

正确选择的参数

让我们尝试将样本的"右"半部分的点拟合成指数概率密度函数。:

有时这种分布是使用单参数函数构造的，但在我们的例子中，如果 .

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]:

第二误差图看起来更均匀。

-值表明我们考虑了更多的系统效应，并且样本和我们编译的函数之间的误差比样本和分布函数之间的误差更接近正态分布的随机变量。 LogNorm.

结论

我们解决了一个应用程序问题，我们需要在示例中拟合几个函数。 因此，我们构建了一个复杂的分布函数来描述物理实验的结果。