广义帕累托分布

广义Pareto (GP)是一个右偏分布，参数化的形状参数，k，和尺度参数，sigma。K也被称为“尾部指数”参数，可以是正的、零的或负的。

x = linspace (0, 1000);情节(x, gppdf (x,。4,1),“- - -”x, x, gppdf (0, 1),“- - -”, x, gppdf (x 2 1),“- - -”）;包含(“x /σ”）;ylabel (的概率密度）;传奇({“k < 0”“k = 0”“k > 0”}）;

注意，当k < 0时，GP超过上限-(1/k)的概率为零。对于k >= 0, GP没有上限。此外，GP通常与第三个阈值参数一起使用，该参数将下限从零移开。这里我们不需要那种一般性。

GP分布是指数分布(k = 0)和帕累托分布(k > 0)的推广。GP将这两种分布包含在一个更大的家族中，因此形状的连续范围是可能的。

模拟超过数数据

GP的分布可以用超越量来建设性地定义。从右尾为零的概率分布开始，比如正态分布，我们可以从该分布中独立抽取随机值。如果我们修复一个阈值，抛出所有低于阈值的值，并从没有抛出的值中减去阈值，结果称为超过。超出部分的分布近似于全科医生。类似地，我们可以在分布的左尾设置一个阈值，并忽略所有高于该阈值的值。阈值必须在原始分布的尾部足够远，以使近似是合理的。

原始分布决定了生成的GP分布的形状参数k。尾部脱落为多项式的分布，如Student的t，会导致一个正的形状参数。尾部呈指数减少的分布，如正态分布，对应的形状参数为零。具有有限尾部的分布，如beta，对应于一个负形状参数。

GP分布在现实世界中的应用包括对股票市场回报的极端情况建模，以及对极端洪水建模。对于本例，我们将使用模拟的数据，这些数据是由5个自由度的Student的t分布生成的。我们将从t分布中提取2000个观察值中最大的5%，然后减去95%分位数以得到超过。

rng (3“旋风”）;x = trnd (2000 1);q =分位数(x, .95);Y = x(x>q) - q;n =元素个数(y)

n = 100

使用最大似然拟合分布

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

paramEsts = gpfit (y);阿拉伯茶= paramEsts (1)尾指数参数sigmaHat = paramEsts (2)%尺度参数

kHat = 0.0987 sigmaHat = 0.7156

正如预期的那样，由于模拟数据是使用t分布生成的，所以k的估计是正的。

视觉检查Fit

为了直观地评估拟合的好坏，我们将绘制一个尾部数据的比例直方图，覆盖我们估计的GP的密度函数。直方图被缩放，使条形图的高度乘以它们的宽度总和为1。

垃圾箱= 0:.25:7;h =酒吧(垃圾箱,histc (y,垃圾箱)/ ((y) * .25长度),“histc”）;h.FaceColor =[。9。9。9);ygrid = linspace(0, 1.1 *马克斯(y), 100);线(ygrid gppdf(阿拉伯茶,ygrid sigmaHat));xlim ([0, 6]);包含(“超过数”）;ylabel (的概率密度）;

我们使用了一个相当小的箱宽，因此在直方图中有大量的噪声。即便如此，拟合的密度仍然遵循数据的形状，所以GP模型似乎是一个不错的选择。

我们也可以比较经验的CDF和拟合的CDF。

[F,彝族]= ecdf (y);情节(咦,gpcdf(咦,阿拉伯茶,sigmaHat),“- - -”）;持有在；楼梯(咦,F,“r”）;持有从；传奇(“拟合广义Pareto CDF”，“经验提供”，“位置”，“东南”）;

参数估计的计算标准误差

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

[nll,acov] = gpllike (paramEsts, y);stdErr =√诊断接头(acov))

stdErr = 0.1158 0.1093

这些标准误差表明，对k的估计的相对精度比对sigma的估计低很多——它的标准误差是估计本身的顺序。形状参数通常很难估计。重要的是要记住，这些标准误差的计算假设GP模型是正确的，并且我们有足够的数据来保持协方差矩阵的渐近逼近。

检验渐近正态性假设

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

然而，这个正态近似可能是，也可能不是一个好的近似。为了评估它在本例中的性能，我们可以使用引导模拟。我们将通过从数据中重新抽样生成1000个重复数据集，为每个数据拟合GP分布，并保存所有重复估计。

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

作为对参数估计量抽样分布的粗略检查，我们可以查看bootstrap重复的直方图。

次要情节(2,1,1);嘘(replEsts (: 1));标题(' k的Bootstrap估计'）;次要情节(2,1,2);嘘(replEsts (:, 2));标题(“sigma的Bootstrap估计”）;

使用参数转换

对于k的bootstrap估计的直方图似乎只有一点不对称，而对于sigma估计的直方图显然是向右倾斜的。对于这种偏度的常见补救方法是在对数尺度上估计参数及其标准误差，在对数尺度上，正态近似可能更合理。Q-Q图是比直方图更好的评估正态性的方法，因为非正态性显示为不是近似沿着直线的点。我们来检查一下log变换是否合适。

次要情节(1、2、1);qqplot (replEsts (: 1));标题(' k的Bootstrap估计'）;次要情节(1、2、2);qqplot(日志(replEsts (:, 2)));标题(“log(sigma)的Bootstrap估计”）;

对k和log(sigma)的自举估计似乎可以接受地接近正态。在未记录的尺度上，对sigma的估计的Q-Q图将证实我们已经在直方图中看到的偏度。因此，更合理的做法是，在正态假设下，首先计算log(sigma)的置信区间为1，然后取指数将该区间转换为sigma的原始尺度。

实际上，这个函数就是这样的gpfit在幕后。

[paramEsts, paramCI] = gpfit (y);

kHat kCI = paramCI(:，1)

kHat = 0.0987 kCI = -0.1283 0.3258

sigmaHat = paramCI(:，2)

sigmaHat = 0.7156 sigmaCI = 0.5305 0.9654

注意，虽然k的95%置信区间是关于最大似然估计对称的，但sigma的置信区间不是。这是因为它是通过变换log()的对称CI而产生的。