广义帕累托分布

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

X = linspace(0,10,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's t分布生成。我们将从t分布中取出2000个观测值中最大的5%，然后减去95%分位数以得到超出值。

rng (3“旋风”）;X = trnd(5,2000,1);Q =分位数(x，.95);Y = x(x>q) - q;N =数字(y)

N = 100

利用最大似然拟合分布

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

参数= gpfit(y);kHat =参数值(1)尾指数参数sigmaHat =参数(2)比例参数

kHat = 0.0987 sigmaHat = 0.7156

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

目测是否合身

为了直观地评估拟合有多好，我们将绘制尾部数据的缩放直方图，与我们估计的GP的密度函数叠加。直方图按比例缩放，使条形图的高度乘以宽度之和为1。

bin = 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);线(ygrid gppdf(阿拉伯茶,ygrid sigmaHat));xlim ([0, 6]);包含(“超过数”）;ylabel (的概率密度）;

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

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

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

计算参数估计的标准误差

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

[nll,acov] = gplike(参数值，y);stdErr = sqrt(diag(acov))

stdErr = 0.1158 0.1093

这些标准误差表明，k估计的相对精度比西格玛估计的相对精度要低一些——其标准误差是估计本身的数量级。形状参数通常难以估计。重要的是要记住，这些标准误差的计算假设GP模型是正确的，并且我们有足够的数据来进行协方差矩阵的渐近逼近。

检验渐近正态性假设

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

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

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

作为参数估计器抽样分布的粗略检查，我们可以查看自举复制的直方图。

次要情节(2,1,1);嘘(replEsts (: 1));标题(k的自举估计）;次要情节(2,1,2);嘘(replEsts (:, 2));标题(“西格玛的引导估计”）;

使用参数转换

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

次要情节(1、2、1);qqplot (replEsts (: 1));标题(k的自举估计）;次要情节(1、2、2);qqplot(日志(replEsts (:, 2)));标题(“对数(西格玛)的引导估计”）;

k和log(sigma)的自举估计似乎可以接受地接近正态。在未记录的尺度上，sigma估计的Q-Q图将确认我们已经在直方图中看到的偏态。因此，在正态性的假设下，首先计算log(sigma)的置信区间，然后对该区间进行幂运算，将其转换回sigma的原始尺度，从而构建σ的置信区间将更加合理。

事实上，这就是这个函数gpfit做幕后的事。

[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变换而产生的。