原文链接:http://tecdat.cn/?p=22566 

原文出处:拓端数据部落公众号

 

本文是极端值推断的内容。我们在广义帕累托分布上使用最大似然方法。

  • 极大似然估计

在参数模型的背景下,标准技术是考虑似然的最大值(或对数似然)。考虑到一些技术性假设,如 R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用 ,R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_02的某个邻域,那么

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_03

其中R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_04表示费雪信息矩阵。在此考虑一些样本,来自广义帕累托分布,参数为R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_05 ,因此 

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_06

 

如果我们解决极大似然的一阶条件,我们得到一个满足以下条件的估计R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_07

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_08

这种渐进正态性的概念如下:如果样本的真实分布是一个具有参数R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_09的GPD,那么,如果n足够大,就会有一个联合正态分布R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_10。因此,如果我们产生大量的样本(足够大,例如200个观测值),那么估计的散点图应该与高斯分布的散点图相同。

  1.  
     
  2.  
    > for(s in 1:1000){
  3.  
    + param[s,]=gpd(x,0)$par.ests
  4.  
     
  5.  
     
  6.  
    > image(x,y,z)

得到一个3D的表示

  1.  
     
  2.  
    > persp(x,y,t(z)
  3.  
    + xlab="xi",ylab="sigma")

 

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_11 R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_12

 

有了200个观测值,如果真正的基础分布是GPD,那么,联合分布R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_13是正态的。
 

  • Delta德尔塔法

另一个重要的属性是德尔塔法。这个想法是,如果是渐进正态,足够平滑,那么也是渐进高斯的。

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_14

从这个属性中,我们可以得到R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_15(这是极值模型中使用的另一个参数化)的正态性,或者在任何四分位数R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_16上 。我们运行一些模拟,再一次检查联合正态性。

  1.  
     
  2.  
    > for(s in 1:1000)
  3.  
    + gpd(x,0)$par.ests
  4.  
    + q=sha * (.01^(-xih) - 1)/xih
  5.  
    + tvar=q+(sha + xih * q)/(1 - xih)
  6.  
    dmnorm(cbind(vx,vy),m,S)
  7.  
    > image(x,y,t(z)

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_17

正如我们所看到的,在样本大小为200的情况下,我们不能使用这个渐进式的结果:看起来我们没有足够的数据。但是,如果我们在n=5000运行同样的代码,

 
  1.  
    > n=5000
  2.  
     
R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_18  

我们得到R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_19R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_20的联合正态性。这就是我们可以从这个结果中得到的delta-方法。

 

  • 轮廓似然( Profile Likelihood )

另一个有趣的方法是Profile 似然函数的概念。因为尾部指数R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_21R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_22在这里是辅助参数。
这可以用来推导出置信区间。在GPD的情况下,对于每个R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_23 ,我们必须找到一个最优的 。我们计算Profile 似然函数,即R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_24 。而我们可以计算出这个轮廓似然的最大值。一般来说,这个两阶段的优化与(全局)最大似然是不等价的,计算结果如下

 

  1.  
     
  2.  
    + profilelikelihood=function(beta){
  3.  
    + -loglik(XI,beta) }
  4.  
    + L[i]=-optim(par=1,fn=profilelik)$value }
  5.  
     

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_25

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_26

 

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_27

 

如果我们想计算轮廓似然的最大值(而不是像以前那样只计算网格上的轮廓似然的值),我们使用

  1.  
     
  2.  
    + profile=function(beta){
  3.  
    + -loglikelihood(XI,beta) }
  4.  
    (OPT=optimize(f=PL,interval=c(0,3)))
  5.  
     

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_28

我们得到结果和最大似然估计的R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_29相似。我们可以用这种方法来计算置信区间,在图表上将其可视化

  1.  
     
  2.  
    > line(h=-up-qchisq(p=.95,df=1)
  3.  
    > I=which(L>=-up-qchisq(p=.95,df=1))
  4.  
    > lines(XIV[I]

R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_30

 

竖线是参数R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_3195%置信区间的下限和上限。


R语言极值推断:广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法_R语言极值推断:广义帕累托分布GPD使用_32

最受欢迎的见解

1.R语言POT超阈值模型和极值理论分析

2.R语言极值理论EVT:基于GPD模型的火灾损失分布分析

3.R语言有极值(EVT)依赖结构的马尔可夫链(MC)对洪水极值分析

4.R语言回归中的hosmer-lemeshow拟合优度检验

5.matlab实现MCMC的马尔可夫切换ARMA – GARCH模型估计

6.R语言区间数据回归分析

7.R语言WALD检验 VS 似然比检验

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI,NRI指标