Engee 文档
Notebook

通过部分抽样确定分布参数

在本例中，我们将把样本分成两部分，并用不同的概率密度函数对每一半进行近似，使用最小二乘法求出每个函数的参数。

创建样本

我们的样本将是线段 [-20,20] 上的布朗过程模型。粒子向左和向右过渡的概率相同，运行长度由给定单位均值的指数分布决定。

让我们从该区域的中心出发，发射 100000 条轨迹，得到一组该过程的现实。

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(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 之间的误差更接近于正态分布的随机变量。

结论

我们解决了一个需要将多个函数拟合到样本中的应用问题。因此，我们构建了一个复杂的分布函数来描述物理实验的结果。