立即登录 马上注册
首页 › 开发帮助文档 › 编程语言 › R语言极值推断：广义帕累托分布GPD使用极大似然估计、轮 ...

赞

评论

收藏

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

编程语言 编程语言 发布于：2021-07-01 21:11 | 阅读数：603 | 评论：0

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

原文出处：拓端数据部落公众号
　　本文是极端值推断的内容。我们在广义帕累托分布上使用最大似然方法。

  • 极大似然估计
　　在参数模型的背景下，标准技术是考虑似然的最大值（或对数似然）。考虑到一些技术性假设，如  DSC0000.png  ， DSC0001.png 的某个邻域，那么
DSC0002.png

　　其中 DSC0003.png 表示费雪信息矩阵。在此考虑一些样本，来自广义帕累托分布，参数为 DSC0004.png  ，因此 
DSC0005.png

　　如果我们解决极大似然的一阶条件，我们得到一个满足以下条件的估计 DSC0006.png
DSC0007.png

　　这种渐进正态性的概念如下：如果样本的真实分布是一个具有参数 DSC0008.png 的GPD，那么，如果n足够大，就会有一个联合正态分布 DSC0009.png 。因此，如果我们产生大量的样本（足够大，例如200个观测值），那么估计的散点图应该与高斯分布的散点图相同。
[list=1]
[*]
　　
[*]
　　> for(s in 1:1000){
　　
[*]
　　+ param[s,]=gpd(x,0)$par.ests
　　
[*]
　　
[*]
　　
[*]
　　> image(x,y,z)
　　
[/list]
　　得到一个3D的表示
[list=1]
[*]
　　
[*]
　　> persp(x,y,t(z)
　　
[*]
　　+ xlab="xi",ylab="sigma")
　　
[/list]
　　
DSC00010.png DSC00011.png

　　有了200个观测值，如果真正的基础分布是GPD，那么，联合分布 DSC00012.png 是正态的。
 

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

　　从这个属性中，我们可以得到 DSC00014.png （这是极值模型中使用的另一个参数化）的正态性，或者在任何四分位数 DSC00015.png 上 。我们运行一些模拟，再一次检查联合正态性。
[list=1]
[*]　　
[*]
　　> for(s in 1:1000)
　　
[*]
　　+ gpd(x,0)$par.ests
　　
[*]
　　+ q=sha * (.01^(-xih) - 1)/xih
　　
[*]
　　+ tvar=q+(sha + xih * q)/(1 - xih)
　　
[*]
　　dmnorm(cbind(vx,vy),m,S)
　　
[*]
　　> image(x,y,t(z)
　　
[/list]
DSC00016.png

　　正如我们所看到的，在样本大小为200的情况下，我们不能使用这个渐进式的结果：看起来我们没有足够的数据。但是，如果我们在n=5000运行同样的代码， 
[list=1]
[*]
　　> n=5000
　　
[*]
　　
[/list]
　　
DSC00017.png  

　　我们得到 DSC00018.png 和 DSC00019.png 的联合正态性。这就是我们可以从这个结果中得到的delta-方法。

 

  • 轮廓似然( Profile Likelihood )
　　另一个有趣的方法是Profile 似然函数的概念。因为尾部指数 DSC00020.png ， DSC00021.png 在这里是辅助参数。
这可以用来推导出置信区间。在GPD的情况下，对于每个 DSC00022.png  ，我们必须找到一个最优的 。我们计算Profile 似然函数，即 DSC00023.png  。而我们可以计算出这个轮廓似然的最大值。一般来说，这个两阶段的优化与（全局）最大似然是不等价的，计算结果如下
[list=1]
[*]　　
[*]
　　+ profilelikelihood=function(beta){
　　
[*]
　　+ -loglik(XI,beta) }
　　
[*]
　　+ L[i]=-optim(par=1,fn=profilelik)$value }
　　
[*]
　　
[/list]
DSC00024.png

DSC00025.png

DSC00026.png

　　如果我们想计算轮廓似然的最大值（而不是像以前那样只计算网格上的轮廓似然的值），我们使用
[list=1]
[*]
　　
[*]
　　+ profile=function(beta){
　　
[*]
　　+ -loglikelihood(XI,beta) }
　　
[*]
　　(OPT=optimize(f=PL,interval=c(0,3)))
　　
[*]
　　
[/list]
DSC00027.png

　　我们得到结果和最大似然估计的 DSC00028.png 相似。我们可以用这种方法来计算置信区间，在图表上将其可视化
[list=1]
[*]
　　
[*]
　　> line(h=-up-qchisq(p=.95,df=1)
　　
[*]
　　> I=which(L>=-up-qchisq(p=.95,df=1))
　　
[*]
　　> lines(XIV[I]
　　
[/list]
DSC00029.png

　　竖线是参数 DSC00030.png 95%置信区间的下限和上限。
DSC00031.jpg

　　最受欢迎的见解
　　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指标

　　
关注下面的标签，发现更多相似文章
免责声明：

1. 本站所有资源来自网络搜集或用户上传，仅作为参考不担保其准确性！

2. 本站内容仅供学习和交流使用，版权归原作者所有！© 查看更多

3. 如有内容侵害到您，请联系我们尽快删除，邮箱：kf@codeae.com

分享到： QQ好友和群QQ好友和群 QQ空间QQ空间

举报

返回列表
相关帖子
avatar
284 经验
1094 文档
OPNE在线字典
  • 手机版| 免责声明|
2021-2024  版权所有：CodeAE代码之家   
粤ICP备2023096576号-2