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

打开脚本

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

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

用最小二乘拟合指数分布

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

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

rng(37岁,“旋风”）;N = 100;X = extend (2,n,1);

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

X = sort(X);P = ((1:n)-0.5)' ./ n;楼梯(x, p,“k -”）;包含(“x”）;ylabel (累积概率(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 = y \ x

muHat = 1.8627

绘制数据和拟合直线。

情节(x, y,“+”y * muHat y“r——”）;包含(“x”）;ylabel ('y = -log(1-p)'）;

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

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

muMLE = expfit(x)

muMLE = 1.7894

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

楼梯(x, p,“k -”）;持有在Xgrid = linspace(0,1.1*max(x)，100)';情节(xgrid expcdf (xgrid muHat),“r——”、xgrid expcdf (xgrid muMLE),“b——”）;持有从包含(“x”）;ylabel (累积概率(p)）;传奇({“数据”，“LS配合”，毫升装的}，“位置”，“东南”）;

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

拟合威布尔分布

对于一个稍微复杂一点的例子，模拟威布尔分布的一些样本数据，并计算x的ECDF。

N = 100;X = wblrnd(2,1,n,1);X = sort(X);P = ((1:n)-0.5)' ./ n;

为了拟合这些数据的威布尔分布，请注意威布尔分布的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 = 2.1420 1.0843

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

情节(计算lnx,呆呆的,“+”，log (paramHat(1)) + logy/paramHat(2)，logy，“r——”）;包含(日志(x)的）;ylabel (日志(日志(1 - p))”）;

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

paramMLE = wblfit(x) stairs(x,p，“k”）;持有在Xgrid = linspace(0,1.1*max(x)，100)';情节(xgrid wblcdf (xgrid paramHat (1) paramHat (2)),“r——”，.．.xgrid, wblcdf (xgrid paramMLE (1) paramMLE (2)),“b——”）;持有从包含(“x”）;ylabel (累积概率(p)）;传奇({“数据”，“LS配合”，毫升装的}，“位置”，“东南”）;

paramMLE = 2.1685 1.0372

A阈值参数示例

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

为了说明这一点，我们将模拟带有阈值的三参数威布尔分布的一些数据。如上所述，为了示例的目的，我们假设数据不知道来自特定的模型，阈值也不知道。

N = 100;X = wblrnd(4,2,n,1) + 4;嘘(x, 20);xlim ([0 16]);

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

X = sort(X);P = ((1:n)-0.5)' ./ n;Logy = log(-log(1-p));Logxm1 = log(x-1);Poly1 = polyfit(log(-log(1-p))，log(x-1)，1);paramHat1 = [exp(poly1(2)) 1/poly1(1)] plot(logxm1,logy，“b +”，log (paramHat1(1)) + logy/paramHat1(2)，logy，“r——”）;包含(“日志(x - 1)”）;ylabel (日志(日志(1 - p))”）;

paramHat1 = 7.4305 4.5574

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

Logxm2 = log(x-2);Poly2 = polyfit(log(-log(1-p))，log(x-2)，1);paramHat2 = [exp(poly2(2)) 1/poly2(1)]

paramHat2 = 6.4046 3.7690

Logxm4 = log(x-4);Poly4 = polyfit(log(-log(1-p))，log(x-4)，1);paramHat4 = [exp(poly4(2)) 1/poly4(1)]

paramHat4 = 4.3530 1.9130

情节(logxm1,呆呆的,“b +”logxm2呆呆的,' r + 'logxm4呆呆的,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——”）;包含(log(x - c)）;ylabel ('log(-log(1 - p))'）;传奇({'Threshold = 1''Threshold = 2''阈值= 4'}，“位置”，“西北”）;

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

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

r2 = @ (x, y) 1 -范数(y - polyval (polyfit (x, y, 1), x))。^2 /范数(y -均值(y)).^2;threshObj = @(c) -r2(log(-log(1-p))，log(x-c));cHat = fminbnd(threshObj，.75*min(x)， .9999*min(x));poly = polyfit(log(-log(1-p))，log(x-cHat)，1);paramHat = [exp(poly(2)) 1/poly(1) cHat] logx = log(x-cHat);Logy = log(-log(1-p));情节(计算lnx,呆呆的,“b +”，log (paramHat(1)) + logy/paramHat(2)，logy，“r——”）;包含('log(x - cHat)'）;ylabel ('log(-log(1 - p))'）;

paramHat = 4.7448 2.3839 3.6029 .使用本文

Non-Location-Scale家庭

指数分布是一个尺度族，对数尺度上，威布尔分布是一个位置尺度族，所以最小二乘法在这两种情况下都很简单。拟合位置-规模分布的一般程序是

计算观测数据的ECDF。
变换分布的CDF，得到数据的某个函数和累积概率的某个函数之间的线性关系。这两个函数不涉及分布参数，但涉及直线的斜率和截距。
将ECDF中的x和p值代入转换后的CDF中，并使用最小二乘拟合一条直线。
用直线的斜率和截距来求解分布参数。

我们还看到，拟合一个具有额外阈值参数的位置规模分布家族只是稍微困难一些。

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

如果将来自ECDF的经验概率与来自参数模型的拟合概率绘制在一起，则沿着从0到1的1:1线的紧密分散表明参数值定义了一个很好地解释观测数据的分布，因为拟合的CDF很好地近似于经验CDF。其思想是找到使概率图尽可能接近1:1线的参数值。如果分布不是数据的良好模型，这甚至可能是不可能的。如果P-P图显示系统地偏离1:1的直线，那么模型可能是有问题的。然而，重要的是要记住，由于这些图中的点不是独立的，解释与回归残差图并不完全相同。

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

N = 100;X = gamnd (2,1,n,1);

计算x的ECDF。

X = sort(X);pEmp = ((1:n)-0.5)' ./ n;

我们可以使用gamma分布参数的任何初始猜测来绘制概率图，比如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);Plot ([0 1]，[0 1]，“k——”pEmp p0Fit,“b +”pEmp p1Fit,' r + '）;包含(“经验概率”）;ylabel (“(临时)拟合伽马概率”）;传奇({“1:1线”，“a = 1, b = 1”，“= 2,b = 1”}，“位置”，“东南”）;

如果你用P-P图的直度来衡量“兼容性”，那么a和b的第二组值就能得到更好的图，因此与数据更兼容。

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

wgt = 1 ./√(pEmp.*(1-pEmp));gammaObj = @ (params)和(重量。* (gamcdf (x, exp (params (1)), exp (params (2))) -pEmp)。^ 2);paramHat = fminsearch(gammaObj，[log(a1)，log(b1)]);paramHat = exp(paramHat)

paramHat = 2.2759 0.9059

pFit = gamcdf(x,paramHat(1)，paramHat(2));Plot ([0 1]，[0 1]，“k——”pEmp pFit,“b +”）;包含(“经验概率”）;ylabel (“拟合伽马概率”）;

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

模型Misspecification

P-P图也可用于比较来自不同分布族的拟合。如果我们试图将这些数据拟合成对数正态分布会发生什么?

wgt = 1 ./√(pEmp.*(1-pEmp));LNobj = @ (params)和(重量。* (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 = 0.5331 0.7038

pFitLN = logcdf (x,paramHatLN(1)，paramHatLN(2));持有在情节(pEmp pFitLN,“处方”）;持有从ylabel (“拟合概率”）;传奇({“1:1线”，“安装伽马”，符合对数正态的}，“位置”，“东南”）;

请注意，对数正态拟合与尾部的伽马拟合是如何系统地不同的。左尾的生长速度较慢，右尾的死亡速度较慢。gamma似乎更符合数据。

A对数正态阈值参数示例

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

为了说明这一点，我们将模拟带有阈值的三参数对数正态分布中的一些数据。

N = 200;X = lognrnd(0，.5,n,1) + 10;嘘(x, 20);xlim (15 [8]);

计算x的ECDF，并找到最佳拟合三参数对数正态分布的参数。

X = sort(X);pEmp = ((1:n)-0.5)' ./ n;wgt = 1 ./√(pEmp.*(1-pEmp));LN3obj = @ (params)和(重量。* (logncdf (x-params(3),参数个数(1),exp (params (2))) -pEmp)。^ 2);C0 = .99*min(x);Mu0 = mean(log(x-c0));Sigma0 = std(log(x-c0));paramHat = fminsearch(LN3obj，[mu0,log(sigma0)，c0]);paramHat(2) = exp(paramHat(2))

paramHat = -0.0698 0.5930 10.1045

pFit = logcdf (x-paramHat(3)，paramHat(1)，paramHat(2));情节(pEmp pFit,“b +”，[0 1]，[0 1]，“k——”）;包含(“经验概率”）;ylabel (拟合3参数对数正态概率）;

精度测量

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

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

estsSim = 0 (1000,3);为i = 1:size(estsSim,1) xSim = lognnd (paramHat(1)，paramHat(2)，n,1) + paramHat(3);estsSim(i，:) = logn3fit(xSim);结束性病(estsSim)

Ans = 0.1542 0.0908 0.1303

查看估计的分布，检查近似正态性的假设对于这个样本量是否合理，或者检查偏差，也可能是有用的。

次要情节(3,1,1),嘘(estsSim (: 1), 20);标题(“日志定位参数引导估计”）;次要情节(3、1、2),嘘(estsSim (:, 2), 20);标题(“对数尺度参数引导估计”）;次要情节(3、1,3),嘘(estsSim (:, 3), 20);标题(“阈值参数引导估计”）;

显然，阈值参数的估计量是倾斜的。这是意料之中的，因为它以最小数据值为界。另外两个直方图表明，对于日志位置参数(第一个直方图)，近似正态性可能也是一个有问题的假设。上面计算的标准误差必须考虑到这一点，通常的置信区间结构可能不适用于日志定位和阈值参数。

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

[paramHat;意思是(estsSim)]

Ans = -0.0698 0.5930 10.1045 -0.0690 0.5926 10.0905

最后，我们也可以用这个函数bootstrp计算引导标准误差估计。它们没有对数据做任何参数假设。

estsBoot = bootstrp(1000，@logn3fit,x);性病(estsBoot)

Ans = 0.1490 0.0785 0.1180

自举标准误差与蒙特卡罗计算相差不远。这并不奇怪，因为拟合的模型与生成示例数据的模型是相同的。

总结

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

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