广义帕累托分布

广义帕累托 (GP) 分布是一种右偏态分布，使用形状参数 k 和尺度参数 sigma 进行参数化。k 也称为“尾部指数”参数，可以为正值、零或负值。

x = linspace(0,10,1000);
plot(x,gppdf(x,-.4,1),'-', x,gppdf(x,0,1),'-', x,gppdf(x,2,1),'-');
xlabel('x / sigma');
ylabel('Probability density');
legend({'k < 0' 'k = 0' 'k > 0'});

请注意，当 k < 0 时，GP 高于上限 -(1/k) 的概率为零。当 k >= 0 时，GP 没有上限。此外，GP 还常常与第三个参数即阈值参数结合使用，该阈值参数使下限偏离零。我们这里不需要用到该参数。

GP 分布是指数分布 (k = 0) 和帕累托分布 (k>0) 的广义化。GP 将这两个分布包括在更大的族中，因此可以实现连续的形状范围。

仿真超阈值数据

GP 分布可从超阈值的角度来定义。基于右尾下降到零的概率分布（例如正态分布），我们可以从该分布独立地抽取随机值。如果我们固定一个阈值，排除所有低于该阈值的值，并将未排除的值减去阈值，那么得出的结果就是超出值。超出值的分布近似于 GP。同样，我们可以在分布的左尾设置一个阈值，并忽略高于该阈值的所有值。阈值必须足够远离原始分布的尾部，这样，逼近才会合理。

原始分布决定生成的 GP 分布的形状参数 k。尾部下降为多项式的分布（例如 Student t 分布）得到的形状参数为正值。尾部呈指数级减小的分布（例如正态分布）对应的形状参数为零。具有有限尾部的分布（例如 beta 分布）对应的形状参数为负值。

GP 分布的实际应用包括对股市收益极值进行建模以及对特大洪水进行建模。对于本示例，我们将使用从具有 5 个自由度的 Student t 分布生成的仿真数据。我们将从 t 分布的 2000 个观测值中提取最大的 5% 的值，然后减去 95% 分位数以获得超阈值。

rng(3,'twister');
x = trnd(5,2000,1);
q = quantile(x,.95);
y = x(x>q) - q;
n = numel(y)

使用最大似然进行分布拟合

GP 分布定义为 0 < sigma，-Inf < k < Inf。然而，当 k < -1/2 时，最大似然估计结果的解释存在问题。幸运的是，这些情况对应于 beta 或三角等分布的尾部拟合，所以这里不会出现问题。

paramEsts = gpfit(y);
kHat      = paramEsts(1)   % Tail index parameter
sigmaHat  = paramEsts(2)   % Scale parameter
kHat =
    0.0987
sigmaHat =
    0.7156

正如所预料的那样，由于仿真数据是使用 t 分布生成的，因此 k 的估计值为正。

通过图形检查拟合效果

为了直观地评估拟合效果，我们将绘制尾部数据的缩放直方图，然后叠加我们估计的 GP 的密度函数。因为是缩放直方图，所以条形高度乘以宽度的总和为 1。

bins = 0:.25:7;
h = bar(bins,histc(y,bins)/(length(y)*.25),'histc');
h.FaceColor = [.9 .9 .9];
ygrid = linspace(0,1.1*max(y),100);
line(ygrid,gppdf(ygrid,kHat,sigmaHat));
xlim([0,6]);
xlabel('Exceedance');
ylabel('Probability Density');

我们使用了比较小的 bin 宽度，因此直方图中存在大量噪声。即便如此，拟合后的密度也符合数据的形状，因此 GP 模型看起来是不错的选择。

我们还可以将经验 CDF 与拟合的 CDF 进行比较。

[F,yi] = ecdf(y);
plot(yi,gpcdf(yi,kHat,sigmaHat),'-');
hold on;
stairs(yi,F,'r');
hold off;
legend('Fitted Generalized Pareto CDF','Empirical CDF','location','southeast');

计算参数估计值的标准误差

为了量化估计值的精度，我们将使用从最大似然估计量的渐近协方差矩阵计算的标准误差。函数 gplike 计算协方差矩阵的数值近似值（作为其第二个输出）。我们也可以使用两个输出实参调用 gpfit ，它应该返回形参的置信区间。

[nll,acov] = gplike(paramEsts, y);
stdErr = sqrt(diag(acov))
stdErr =
    0.1158
    0.1093

这些标准误差表明 k 估计值的相对精度远低于 sigma 的相对精度（sigma 的标准误差与估计值本身相似）。形状参数通常难以估计。必须要记住的是，计算这些标准误差的前提是假设 GP 模型正确，并且我们有足够的数据保证对协方差矩阵的渐近逼近是合理的。

检查渐近正态性假设

对标准误差的解释通常基于如下假设：如果可以对来自相同数据源的数据多次重复相同的拟合，那么参数的最大似然估计值将大致遵循正态分布。例如，置信区间通常基于这一假设。

然而，正态逼近可能好，也可能不好。为了评估它在此示例中的良好程度，我们可以使用自助仿真。我们通过重抽样从数据中生成 1000 个副本数据集，然后对每个数据集进行 GP 分布拟合，并保存所有副本估计值。

replEsts = bootstrp(1000,@gpfit,y);

要大致检查参数估计量的抽样分布，我们可以查看自助复制数据的直方图。

subplot(2,1,1);
hist(replEsts(:,1));
title('Bootstrap estimates of k');
subplot(2,1,2);
hist(replEsts(:,2));
title('Bootstrap estimates of sigma');

使用参数变换

k 的自助估计值的直方图看起来只有一点不对称，而 sigma 估计值的直方图明显向右偏斜。为了纠正这一偏斜，常见的做法是基于对数尺度来估计参数及其标准误差，因为基于对数尺度的正态逼近可能更为合理。在评估正态性方面，Q-Q 图比直方图更好，因为非正态性的数据会表现为落在一条近似直线之外的点。我们来检验一下，看看对数变换对 sigma 是否合适。

subplot(1,2,1);
qqplot(replEsts(:,1));
title('Bootstrap estimates of k');
subplot(1,2,2);
qqplot(log(replEsts(:,2)));
title('Bootstrap estimates of log(sigma)');

k 和 log(sigma) 的自助估计值显示为大致接近正态。基于非对数尺度的 sigma 估计值 Q-Q 图将印证我们已经在直方图中看到的偏斜。因此，先在假设正态的情况下计算 log(sigma) 的置信区间，然后取幂，将该置信区间变换回 sigma 的原始尺度，这样构造 sigma 的置信区间更为合理。

事实上，这正是函数 gpfit 在幕后所做的工作。

[paramEsts,paramCI] = gpfit(y);
kCI  = paramCI(:,1)
kHat =
    0.0987
kCI =
   -0.1283
    0.3258

sigmaHat
sigmaCI  = paramCI(:,2)
sigmaHat =
    0.7156
sigmaCI =
    0.5305
    0.9654

请注意，虽然 k 的 95% 的置信区间围绕最大似然估计都是对称的，但 sigma 的置信区间并不是这样。那是因为它是通过变换 log(sigma) 的对称 CI 创建的。