用最小二乘拟合指数分布

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

首先，模拟一些样本数据。我们将使用指数分布来生成数据。在本例中，与实际情况一样，我们将假设数据不来自特定的模型。

rng（37，“龙卷风”); n=100；x=exprnd（2，n，1）；

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

x = (x);P = ((1:n)-0.5)' ./ n;楼梯(x, p,“k-”)；xlabel(“x”); 伊拉贝尔(‘累积概率（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 vs.-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——”)；xlabel(“x”); 伊拉贝尔('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 * 100 (max (x)) ';情节(xgrid expcdf (xgrid muHat),“r——”、xgrid expcdf (xgrid muMLE),“b——”）;持有从包含(“x”); 伊拉贝尔(‘累积概率（p）’); 传奇({“数据”,“它合适吗？”,“ML-Fit”}，“位置”,“东南”）;

这两种方法给出了非常相似的拟合分布，尽管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））；多边形=多边形拟合（逻辑，对数x，1）；paramHat=[exp（poly（2））1/poly（1）]

参数=2.1420 1.0843

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

绘图（logx、逻辑、，'+'，log (paramHat(1)) + logy/paramHat(2)，“r——”)；xlabel('日志（x）'); 伊拉贝尔('日志（-log（1-p））'）;

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

paramMLE=wblfit（x）楼梯（x，p，“k”）;持有在…上xgrid=linspace（0,1.1*max（x），100）”；绘图（xgrid，wblcdf（xgrid，paramHat（1），paramHat（2）），“r——”,...xgrid，wblcdf（xgrid，参数（1），参数（2）），“b——”）;持有从包含(“x”); 伊拉贝尔(‘累积概率（p）’); 传奇({“数据”,“它合适吗？”,“ML-Fit”}，“位置”,“东南”）;

paramMLE=2.1685 1.0372

A阈值参数示例

有时需要用阈值参数拟合正分布，如Weibull分布或对数正态分布。例如，威布尔随机变量的值超过（0，Inf），阈值参数c将该范围移到（c，Inf）。如果阈值参数已知，则不存在困难。但如果阈值参数未知，则必须对其进行估计。这些模型很难与最大似然法相匹配——这种似然法可能有多种模式，甚至对于数据不合理的参数值，这种似然法变得无限大，因此最大似然法通常不是一种好方法。但是，在最小二乘法的基础上，我们可以得到稳定的估计。

为了说明这一点，我们将模拟一些来自三参数威布尔分布的数据，其中包含一个阈值。如上所述，我们将在示例中假设数据不来自特定模型，并且阈值未知。

n=100；x=wblrnd（4,2，n，1）+4；hist（x，20）；xlim（[016]）；

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

x=排序（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）；参数1=[exp（poly1（2））1/poly1（1）]图（logxm1，逻辑，“b+”，log (paramHat1(1)) + logy/paramHat1(2)，“r——”)；xlabel(“日志(x - 1)”); 伊拉贝尔('日志（-log（1-p））'）;

参数1=7.4305 4.5574

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

logxm2 =日志(x - 2);poly2 = polyfit(日志(日志(1 - p)),日志(x - 2), 1);paramHat2 = [exp(poly2(2)) 1/poly2(1)]

参数2=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+”,...日志（参数1（1））+逻辑/参数1（2），逻辑，“b——”,...日志（paramHat2（1））+Logic/paramHat2（2），Logic，“r——”,...日志（paramHat4（1））+Logic/paramHat4（2），逻辑，“g——”)；xlabel('日志（x-c）'); 伊拉贝尔('log（-log（1-p））'); 传奇({'阈值=1'“阈值= 2”“阈值=4”}，“位置”,“西北”）;

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-范数（y-polyval）（polyfit（x，y，1），x））。^2/范数（y-均值（y））。^2；阈值j=@（c）-r2（log（-log（1-p）），log（x-c））；cHat=fminbnd（阈值，.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=logy=logy（x-cHat））；logy=logy=logy，logy，“b+”，log (paramHat(1)) + logy/paramHat(2)，“r——”)；xlabel('日志（x-聊天）'); 伊拉贝尔('log（-log（1-p））'）;

参数=4.7448 2.3839 3.6029

Non-Location-Scale家庭

指数分布是一个尺度族，而在对数尺度上，威布尔分布是一个位置尺度族，因此在这两种情况下，这种最小二乘法非常简单。拟合位置尺度分布的一般步骤如下

我们还看到，使用额外的阈值参数拟合位置比例族的分布只稍微困难一些。

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

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

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

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

计算x的ECDF。

x = (x);pEmp = ((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）；绘图（[01]，[01]，“k——”，pEmp，p0Fit，“b+”，pEmp，p1Fit，“r+”)；xlabel(“经验概率”); 伊拉贝尔(“(临时)拟合伽玛概率”); 传奇({“1:1线”,“a=1，b=1”,“a=2，b=1”}，“位置”,“东南”）;

a和b的第二组值有助于绘制更好的曲线图，因此，如果通过P-P曲线图的直线度来测量“兼容”，则与数据更兼容。

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

wgt=1./sqrt（pEmp.*（1-pEmp））；gammaObj=@（params）和（wgt.*（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））；绘图（[01]，[01]，“k——”，pEmp，pFit，“b+”)；xlabel(“经验概率”); 伊拉贝尔(“拟合伽马概率”）;

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

模型Misspecification

P-P图对于比较不同分布族的拟合也很有用。如果我们尝试对这些数据拟合对数正态分布，会发生什么？

wgt = 1 ./ sqrt(pmp .*(1- pmp));LNobj = @ (params)和(重量。* (logncdf (x, params (1), exp (params (2))) -pEmp)。^ 2);mu0 =意味着(日志(x));sigma0 =性病(日志(x));paramHatLN = fminsearch (LNobj [mu0、日志(sigma0)]);paramHatLN (2) = exp (paramHatLN (2))

paramHatLN=0.5331 0.7038

pFitLN=logncdf（x，paramHatLN（1），paramHatLN（2））；持有在…上情节(pEmp pFitLN,“处方”）;持有从伊拉贝尔(“拟合概率”); 传奇({“1:1线”,“适配伽马”,符合对数正态的}，“位置”,“东南”）;

请注意对数正态拟合与gamma拟合在尾部的系统性差异。它在左尾部增长较慢，在右尾部消亡较慢。gamma似乎更适合数据。

A日志正常阈值参数示例

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

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

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

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

x=sort（x）；pEmp=（（1:n）-0.5）./n；wgt=1./sqrt（pEmp.*（1-pEmp））；LN3obj=@（params）和（wgt.*（logncdf（x-params（3），params（1），exp（params（2））-pEmp.^2；c0=.99*min（x）；mu0=平均值（log（x-c0））；sigma0=std（log（x-c0））；paramHat=fminsearch（LN3obj，[mu0，log（sigma0），log（sigma0），exp）（param2））

paramHat = -0.0698 0.5930 10.1045

pFit=logncdf（x-paramHat（3）、paramHat（1）、paramHat（2））；绘图（pEmp、pFit、，“b+”, [0 1],[0 1],“k——”)；xlabel(“经验概率”); 伊拉贝尔(“拟合3参数对数正态概率”）;

精密度测量

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

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

estsSim = 0 (1000 3);对于i=1:size（estsSim，1）xSim=lognrnd（paramHat（1），paramHat（2），n，1）+paramHat（3）；estsSim（i，：）=logn3fit（xSim）；终止性病(estsSim)

ans=0.15420.0908 0.1303

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

子地块（3,1,1），历史（ESTSIM（：，1），20）；标题(“日志位置参数引导估计”)第（3,1,2）小批，历史（ESTSIM（：，2），20）；标题(“对数刻度参数引导估计”)第（3,1,3）小批，历史（ESTSIM（：，3），20）；标题(“阈值参数引导估计”）;

显然，阈值参数的估计器是倾斜的。这是意料之中的，因为它是由最小数据值限定的。另外两个柱状图表明，对于对数位置参数（第一个柱状图），近似正态性可能也是一个有问题的假设。在解释上述计算的标准误差时，必须牢记这一点，置信区间的通常构造可能不适合于测井位置和阈值参数。

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

[参数；平均值（estsSim）]

ans=-0.0698 0.5930 10.1045-0.0690 0.5926 10.0905

最后，我们还可以使用这个函数靴筒计算bootstrap标准误差估计。它们不做任何关于数据的参数假设。

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

ans=0.1490.0785 0.1180

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

总结

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

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