Syslab使用手册
⋮
2026a

# 使用累积概率拟合单变量分布

此示例说明如何使用累积分布函数的最小二乘估计来拟合单变量分布。 这是一种普遍适用的方法，在最大似然失败的情况下很有用，例如某些包含阈值参数的模型。

将单变量分布拟合到数据的最常用方法是最大似然法。 但最大似然并非在所有情况下都有效，有时需要使用其他估计方法，例如矩量法。 适用时，最大似然可能是更好的选择，因为它通常更有效。 但是这里描述的方法提供了另一种工具，可以在需要时使用。

# 使用最小二乘法拟合指数分布

“最小二乘法”最常用于拟合回归线或曲面以将响应变量建模为一个或多个预测变量的函数。 这里描述的方法是最小二乘法的一个非常不同的应用：单变量分布拟合，只有一个变量。

首先，首先模拟一些示例数据。 我们将使用指数分布来生成数据。 出于本示例的目的，在实践中，我们假设数据不知道来自特定模型。

using TyPlot
using TyStatistics
using TyMath
using TyOptimization
rng = MT19937ar(37)
n = 100
x = exprnd(rng,2,n)

接下来，计算数据的经验累积分布函数 (ECDF)。 这只是一个阶跃函数，每个数据点 x 的累积概率 p 为 1/n。

figure()
x = sort(x)
p = collect(((1:n).-0.5) / n)
stairs(x, p; color="k")
xlabel("x")
ylabel("Cumulative probability (p)")

我们将为这些数据拟合指数分布。 一种方法是找到其累积分布函数 (CDF) 最接近（在下文中解释的某种意义上）数据的 ECDF 的指数分布。 指数 CDF 是 p = Pr{X <= x} = 1 - exp(-x/mu)。 将其转换为 -log(1-p)*mu = x 可得出 -log(1-p) 和 x 之间的线性关系。 如果数据确实来自指数，如果我们将计算出的 ECDF 中的 x 和 p 值代入该方程，我们应该至少近似地看到线性关系。 如果我们使用最小二乘法来拟合一条通过原点到 x 与 -log(1-p) 的直线，则该拟合线表示“最接近”数据的指数分布。 直线的斜率是参数 mu 的估计值。

同样，我们可以将 y = -log(1-p) 视为来自标准（均值 1）指数分布的“理想化样本”。 这些理想化值在概率标度上完全等距分布。 如果数据来自指数分布，则 x 和 y 的 Q-Q 图应该近似线性，我们将通过原点拟合最小二乘线到 x 与 y。

y =@. -log(1-p)
muHat = reshape(y,100,1) \ reshape(x,100,1)

muHat = 1×1 Matrix{Float64}:
 1.8626678210882612

绘制数据和拟合线。

figure()
plot(x,y,"+", y*muHat,y,"r--")
xlabel("x")
ylabel("y = -log(1-p)")

请注意，我们所做的线性拟合最小化了水平或“x”方向的误差平方和。 这是因为 y = -log(1-p) 的值是确定性的，而 x 值是随机的。 也可以对 y 与 x 进行回归，或使用其他类型的线性拟合，例如加权回归、正交回归，甚至稳健回归。 我们不会在这里探讨这些可能性。

为了进行比较，按最大似然法拟合数据。

muMLE, = expfit(x)

muMLE = 1.7894257093567112

现在在未转换的累积概率尺度上绘制两个估计分布。

figure()
stairs(x,p,color="k")
hold("on")
xgrid = LinRange(0,1.1*maximum(x),100)
plot(xgrid,expcdf.(xgrid,muHat),"r--", xgrid,expcdf(xgrid,muMLE),"b--")
hold("off")
xlabel("x")
ylabel("Cumulative Probability (p)")
legend(["Data","LS Fit","ML Fit","location","southeast"])

这两种方法给出了非常相似的拟合分布，尽管 LS 拟合更多地受到分布尾部观测值的影响。

# 拟合威布尔分布

对于稍微复杂一点的示例，模拟来自 Weibull 分布的一些样本数据，并计算 x 的 ECDF。

n = 100
x = wblrnd(rng,1,2,n)
x = sort(x)
p = collect(((1:n).-0.5)) / n

要将 Weibull 分布拟合到这些数据，请注意 Weibull 的 CDF 是 p = Pr{X <= x} = 1 - exp(-(x/a)^b)。 将其转换为 log(a) + log(-log(1-p))*(1/b) = log(x) 再次给出线性关系，这次是 log(-log(1-p)) 和 log （X）。 我们可以使用最小二乘法来使用 ECDF 中的 p 和 x 在转换后的尺度上拟合一条直线，该直线的斜率和截距会导致对 a 和 b 的估计。

logx = log.(x)
logy =@. log(-log(1 - p))
poly, = polyfit(logy,logx,1)
paramHat = [exp(poly[2]) 1/poly[1]]

paramHat = 1×2 Matrix{Float64}:
 2.14198  1.08431

在变换后的比例上绘制数据和拟合线。

figure()
plot(logx,logy,"+", log(paramHat[1]) .+ logy/paramHat[2],logy,"r--");
xlabel("log(x)")
ylabel("log(-log(1-p))")

为了进行比较，用最大似然法拟合数据，并在未转换的尺度上绘制两个估计分布。

paramMLE, = wblfit(x)
figure()
stairs(x,p,color="k")
hold("on")
xgrid = LinRange(0,1.1*maximum(x),100)'
plot(xgrid,wblcdf(xgrid,paramHat[1],paramHat[2]),"r--",xgrid,wblcdf(xgrid,paramMLE[1],paramMLE[2]),"b--")
hold("off")
xlabel("x")
ylabel("Cumulative Probability (p)")
legend(["Data","LS Fit","ML Fit"],loc="southeast")

paramMLE = 1×2 Matrix{Float64}:
 2.16851  1.03715

# 阈值参数示例

有时需要使用阈值参数来拟合正分布，例如 Weibull 或对数正态分布。 例如，Weibull 随机变量取值超过 (0,Inf)，阈值参数 c 将该范围移动到 (c,Inf)。 如果阈值参数已知，则没有困难。 但如果阈值参数未知，则必须对其进行估计。 这些模型很难用最大似然拟合--似然可以有多种模式，甚至对于数据不合理的参数值变得无穷大，因此最大似然通常不是一个好的方法。 但是通过对最小二乘法进行少量添加，我们可以获得稳定的估计。

为了说明，我们将模拟来自具有阈值的三参数 Weibull 分布的一些数据。 如上所述，为了示例的目的，我们假设数据不知道来自特定模型，并且阈值未知。

n = 100
x = wblrnd(rng,2,4,n,1) .+ 4
figure()
hist(x,20)
xlim([0 16])

我们如何将三参数威布尔分布拟合到这些数据中？如果我们知道阈值是什么，例如1，我们可以从数据中减去该值，然后使用最小二乘法来估计威布尔形状和尺度参数。

x = sort(vec(x))
p = vec(((1:n).-0.5)') / n
logy =@. log(-log(1-p))
logxm1 =@. log(x-1)
poly1, = polyfit(logy,logxm1,1)
paramHat1 = [exp(poly1[2]) 1/poly1[1]]
figure()
plot(logxm1,logy,"b+", log(paramHat1[1]) .+ logy/paramHat1[2],logy,"r--")
xlabel("log(x-1)")
ylabel("log(-log(1-p))")

paramHat1 = 1×2 Matrix{Float64}:
 7.43052  4.55741

这不是很好的拟合--log(x-1) 和 log(-log(1-p)) 没有线性关系。当然，这是因为我们不知道正确的阈值。如果我们尝试减去不同的阈值，我们会得到不同的曲线图和不同的参数估计。

logxm2 =@. log(x-2)
poly2, = polyfit(logy,logxm2,1)
paramHat2 = [exp(poly2[2]) 1/poly2[1]]

paramHat2 = 1×2 Matrix{Float64}:
 6.40463  3.76899

logxm4 =@. log(x-4)
poly4, = polyfit(logy,logxm4,1)
paramHat4 = [exp(poly4[2]) 1/poly4[1]]

paramHat4 = 1×2 Matrix{Float64}:
 4.35297  1.91298

figure()
plot(logxm1,logy,"b+", logxm2,logy,"r+", logxm4,logy,"g+", log(paramHat1[1]) .+ logy/paramHat1[2],logy,"b--", log(paramHat2[1]) .+ logy/paramHat2[2],logy,"r--", log(paramHat4[1]) .+ logy/paramHat4[2],logy,"g--")
xlabel("log(x - c)")
ylabel("log(-log(1 - p))")
legend(["Threshold = 1" "Threshold = 2" "Threshold = 4"], loc="northwest")

log(x-4) 和 log(-log(1-p)) 之间的关系近似线性。 如果我们减去真正的阈值参数，我们希望看到一个近似线性的图，这证明 4 可能是一个合理的阈值。 另一方面，2 和 3 的图更系统地不同于线性图，这证明这些值与数据不一致。

这个论点可以形式化。 对于阈值参数的每个临时值，相应的临时 Weibull 拟合可以表征为使变换变量 log(x-c) 和 log(-log(1-p)) 上的线性回归的 R^2 值最大化的参数值。为了估计阈值参数，我们可以更进一步，并在所有可能的阈值上最大化 R^2 值。

r2(x,y) = 1.0 - norm(y .- polyval(vec(polyfit(x,y,1)[1]),x)).^2 / norm(y .- mean(y))^2
threshObj(c) = -r2(logy,log.(x.-c))
cHat, = fminbnd(threshObj, 0.75 * minimum(x), 0.9999 * minimum(x))
logx =@. log(x-cHat)
logy =@. log(-log(1-p))
poly, = polyfit(logy,logx,1)
paramHat = [exp(poly[2]) 1/poly[1] cHat]

paramHat = 1×3 Matrix{Float64}:
 4.74483  2.38388  3.6029

figure()
plot(logx,logy,"b+", log(paramHat[1]) .+ logy/paramHat[2],logy,"r--")
xlabel("log(x - cHat)")
ylabel("log(-log(1 - p))")

# 非位置尺度族

指数分布是一个尺度族，而在对数尺度上，威布尔分布是一个位置尺度族，所以在这两种情况下，这种最小二乘法是直截了当的。 拟合位置尺度分布的一般过程是

  • 计算观测数据的 ECDF。

  • 变换分布的 CDF 以获得数据的某些函数与累积概率的某些函数之间的线性关系。 这两个函数不涉及分布参数，但直线的斜率和截距涉及。

  • 将 ECDF 中的 x 和 p 值插入转换后的 CDF，并使用最小二乘法拟合一条直线。

  • 根据直线的斜率和截距求解分布参数。

我们还看到，用一个额外的阈值参数来拟合一个位置尺度族的分布只是稍微困难一点。

但其他非位置尺度族的分布，如伽玛分布，有点棘手。 没有 CDF 的转换会给出线性关系。 然而，我们可以使用类似的想法，只是这次在未转换的累积概率尺度上工作。 P-P 图是可视化拟合过程的合适方式。

如果将 ECDF 的经验概率与参数模型的拟合概率作图，则沿着 1:1 线从零到一的紧密散布表明参数值定义的分布很好地解释了观察到的数据，因为拟合 CDF 近似 经验 CDF 很好。 这个想法是找到使概率图尽可能接近 1:1 线的参数值。 如果分布不是数据的良好模型，那甚至可能是不可能的。 如果 P-P 图显示系统偏离 1:1 线，则模型可能有问题。 但是，请务必记住，由于这些图中的点不是独立的，因此解释与回归残差图并不完全相同。

例如，我们将模拟一些数据并拟合伽玛分布。

n = 100
x = gamrnd(rng,2,1,n)

计算 x 的 ECDF 值

x = sort(x)
pEmp = vec(((1:n).-0.5)) / n

我们可以使用伽马分布参数的任何初始猜测来绘制概率图，比如 a=1 和 b=1。 这个猜测不是很好--参数 CDF 的概率与 ECDF 的概率并不接近。 如果我们尝试不同的 a 和 b，我们会在 P-P 图上得到不同的散点，与 1:1 线有不同的差异。 由于我们知道此示例中的真实 a 和 b，因此我们将尝试这些值。

a0 = 1
b0 = 1
p0Fit = gamcdf.(x,a0,b0)
a1 = 2
b1 = 1
p1Fit = gamcdf.(x,a1,b1)
figure()
plot([0 1],[0 1],"k--", pEmp,p0Fit,"b+", pEmp,p1Fit,"r+")
xlabel("Empirical Probabilities")
ylabel("(Provisionally) Fitted Gamma Probabilities")
legend(["1:1 Line","a=1, b=1", "a=2, b=1"], loc="southeast")

a 和 b 的第二组值可以绘制出更好的图，因此与数据更兼容，如果您通过绘制 P-P 图的直度来衡量“兼容”。

为了使散点尽可能接近 1:1 线，我们可以找到 a 和 b 的值，使到 1:1 线的平方距离的加权和最小。 权重是根据经验概率定义的，在图的中心最低，在极端处最高。 这些权重补偿了拟合概率的方差，它在中值附近最高，在尾部最低。 这个加权最小二乘过程定义了 a 和 b 的估计量。

wgt =@. 1 / sqrt(pEmp*(1-pEmp))
gammaObj(params) =  sum(wgt.*(gamcdf.(x,exp(params[1]),exp(params[2])).-pEmp).^2)
paramHat, = fminsearch(gammaObj, [log(a1), log(b1)])
paramHat = exp.(paramHat)

paramHat = 2-element Vector{Float64}:
 2.2758804218649913
 0.9059013402660073

pFit = gamcdf.(x,paramHat[1],paramHat[2])
figure()
plot([0 1],[0 1],"k--", pEmp,pFit,"b+")
xlabel("Empirical Probabilities")
ylabel("Fitted Gamma Probabilities")

请注意，在前面考虑的位置尺度情况下，我们可以用一条直线拟合来拟合分布。 在这里，与阈值参数示例一样，我们必须迭代地找到最合适的参数值。

# 模型规格错误

P-P 图也可用于比较来自不同分布族的拟合。 如果我们尝试对这些数据拟合对数正态分布会怎样？

wgt =@. 1 / sqrt(pEmp*(1-pEmp))
LNobj(params) = sum(wgt.*(logncdf.(x,params[1],exp(params[2])).-pEmp).^2)
mu0 = mean(log.(x))
sigma0 = std(log.(x))
paramHatLN, = fminsearch(LNobj, [mu0,log(sigma0)])
paramHatLN[2] = exp(paramHatLN[2])
paramHatLN

paramHatLN = 2-element Vector{Float64}:
 0.5330920737368392
 0.7038343028942149

pFitLN = logncdf.(x,paramHatLN[1],paramHatLN[2])
hold("on")
plot(pEmp,pFitLN,"rx")
hold("off")
ylabel("Fitted Probabilities")
legend(["1:1 Line", "Fitted Gamma", "Fitted Lognormal"],loc="southeast")

注意到对数正态拟合与伽马拟合在尾部存在系统性地不同。对数正态拟合在前尾处增长地更慢而在后尾处又减少地更慢。这里伽马对数据的拟合效果更好。

# 对数正态阈值参数示例

对数正态分布很容易通过最大似然拟合，因为一旦将对数变换应用于数据，最大似然与拟合正态分布相同。 但有时需要在对数正态模型中估计阈值参数。 这种模型的可能性是无限的，因此最大似然估计不起作用。 然而，最小二乘法提供了一种进行估计的方法。 由于双参数对数正态分布可以对数转换为位置尺度族，我们可以遵循与前面示例中相同的步骤，该示例显示了使用阈值参数拟合 Weibull 分布。 然而，在这里，我们将在累积概率尺度上进行估计，如前一个显示与伽马分布拟合的示例一样。

为了说明，我们将模拟一些来自三参数对数正态分布的数据，并设置一个阈值。

n = 200
x = lognrnd(rng,0,0.5,n,) .+ 10
figure()
hist(x,20)
xlim([8 15])

计算 x 的 ECDF，并找出最适合三参数对数正态分布的参数。

x = sort(x)
pEmp = ((1:n).-0.5) / n
wgt =@. 1 / sqrt(pEmp*(1-pEmp))
LN3obj(params) = sum(wgt.*(logncdf.(x.-params[3],params[1],exp(params[2])).-pEmp).^2)
c0 = 0.99*minimum(x)
mu0 = mean(log.(x.-c0))
sigma0 = std(log.(x.-c0))
paramHat, = fminsearch(LN3obj, [mu0, log(sigma0), c0])
paramHat[2] = exp(paramHat[2])

paramHat = 3-element Vector{Float64}:
 -0.06984083437642054
  0.5930068999222436
 10.104532699540883

pFit = logncdf.(x.-paramHat[3],paramHat[1],paramHat[2])
figure()
plot(pEmp,pFit,"b+", [0 1],[0 1],"k--")
xlabel("Empirical Probabilities")
ylabel("Fitted 3-param Lognormal Probabilities")

# 精度测量

参数估计只是一部分--模型拟合还需要一些衡量估计精确度的指标，通常是标准误差。 对于最大似然，通常的方法是使用信息矩阵和大样本渐近参数来近似估计器在重复采样上的协方差矩阵。 对于这些最小二乘估计器，不存在这样的理论。

然而，蒙特卡洛模拟提供了另一种估计标准误差的方法。 如果我们使用拟合模型生成大量数据集，我们可以用蒙特卡洛标准差来近似估计量的标准误差。 为简单起见，我们首先定义一个拟合函数 logn3fit。

function logn3fit(x)
    x = sort(x)
    n = length(x)
    pEmp = vec(((1:n) .- 0.5) / n)
    wgt = @. 1 / sqrt(pEmp * (1 - pEmp))
    LN3obj(params) = sum(wgt .* (logncdf.(x .- params[3], params[1], exp(params[2])) .- pEmp) .^ 2)
    c0 = 0.95 * minimum(x)
    mu0 = mean(x .- c0)
    sigma0 = std(x .- c0)
    opts = optimset(;MaxIterations=1000, MaxFunctionEvaluations=2000)
    paramEsts, = fminsearch(LN3obj, [mu0,log(sigma0),c0], opts);
    paramEsts[2] = exp(paramEsts[2])
    return paramEsts
end

estsSim = zeros(1000,3);
for i in 1:1000
    xSim = lognrnd(rng,paramHat[1],paramHat[2],n) .+ paramHat[3]
    estsSim[i,:] = logn3fit(xSim)
end
std(estsSim,dims=1)

ans = 1×3 Matrix{Float64}:
 0.154239  0.0908455  0.130284

查看估计值的分布、检查近似正态假设对于该样本量是否合理或检查偏差也可能很有用。

subplot(3,1,1)
hist(estsSim[:,1],20)
title("Log-Location Parameter Bootstrap Estimates")
subplot(3,1,2)
hist(estsSim[:,2],20)
title("Log-Scale Parameter Bootstrap Estimates")
subplot(3,1,3)
hist(estsSim[:,3],20)
title("Threshold Parameter Bootstrap Estimates")
tightlayout()

显然，阈值参数的估计量是有偏差的。 这是意料之中的，因为它受最小数据值限制。 另外两个直方图表明，近似正态性也可能是对数位置参数（第一个直方图）的可疑假设。 上面计算的标准误差必须考虑到这一点来解释，并且置信区间的通常构造可能不适合日志位置和阈值参数。

模拟估计值的均值接近于用于生成模拟数据的参数值，表明该过程在此样本量下大致无偏，至少对于接近估计值的参数值而言是这样。

[paramHat'; mean(estsSim;dims=1)]

ans = 2×3 Matrix{Float64}:
 -0.0698408  0.593007  10.1045
 -0.0690233  0.592567  10.0905

# 总结

这里描述的拟合方法是最大似然的替代方法，当最大似然不能提供有用的参数估计时，该方法可用于拟合单变量分布。一个重要的应用是拟合包含阈值参数的分布，例如三参数对数正态分布。标准误差比最大似然估计更难计算，因为不存在解析近似，但模拟提供了一种可行的选择。

这里用来说明拟合方法的P-P图本身是有用的，作为拟合单变量分布时缺乏拟合的直观指示。