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

由Kaizong Ye，Liao Bao撰写

最近我们被客户要求撰写关于极值推断的研究报告。在参数模型的背景下，标准技术是考虑似然的最大值（或对数似然）。

极大似然估计

考虑到一些技术性假设，如，的某个邻域，那么

简单说来定义是：

假设们有高维似然函数

$L(\theta_1, \theta_2|x)$

那么

$L(\theta_1|x) = L_{\sup \theta_2}(\theta_1, \theta_2|x)$ 叫做对于 $\theta_1$ 的profile likelihood。这个式子表示你给定任意的 $\theta_1$ ,找到在这个值下面（任意改变 $\theta_2$ ）最大的 $L(\theta_1)$ ： $L(\theta_1, \theta_2=?)$

如果是normal distribution, 对于 $\mu$ 和 $\sigma^2$ ，

$L(\mu|x)_{\text{profile}} = L(\mu, \sigma^2 = \hat{\sigma}_{\mu}^2|x) ;\hat{\sigma}_{\mu}^2 = \frac{1}{n} \sum_i^n(x_i - \mu)^2$

注意和其它的几种似然函数微妙的区别。

estimated likelihood (有时人叫它pseudo likelihood。。随便吧）：先计算你的 $\hat{\theta_2}$ 再带入似然函数。例如，你可以把normal distribution, sample variance代入，然后得到对于mean的estimated likelihood。注意，对比第一个式子，这里我们需要用 $\bar{x}$ 替换 $\hat{\sigma}^2$ 中的 $\mu$

$L_\text{estimated} (\mu|x) = L(\mu, \sigma^2 = \hat{\sigma}^2|x);\hat{\sigma}^2=\frac{1}{n} \sum_i^n (x_i - \bar{x})^2$

或者

$L_\text{estimated} (\mu|x) = L(\mu, \sigma^2 = S^2|x);S^2=\frac{1}{n-1} \sum_i^n (x_i - \bar{x})^2$

还有一个叫做integrated likelihood：我们假设 $L(\theta_1, \theta_2)$ 是一个目标函数，同时假设我们的参数有概率分布的，那么就可以积分得到marginal distribution

$L_\text{integrated}(\mu) = \int L(\mu, \sigma^2) f(\sigma^2) d\sigma^2$

这个个是Bayesian 分析中用的，因为prior就是你的参数分布，必须有这个才能积分。

modified profile likelihood：有些人认为 marginal likelihood才是精确的，但是由于难算，就想用profile likelihood的表达式去近似，最后加上高阶修正；或者利用Laplace approximation理解也可以。要达到的目的就是你的modified profile会更加像marginal likelihood,或者真实likelihood. 个人认为，这就是频率派的人在朝着Bayesian派在Likelihood方法靠拢的做法。

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

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

这种渐进正态性的概念如下：如果样本的真实分布是一个具有参数的GPD，那么，如果n足够大，就会有一个联合正态分布。因此，如果我们产生大量的样本（足够大，例如200个观测值），那么估计的散点图应该与高斯分布的散点图相同。

> for(s in 1:1000){
+ param\[s,\]=gpd(x,0)$par.ests


> image(x,y,z)

可下载资源

完整程序、数据和文档（word）

最受欢迎的见解

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指标

得到一个3D的表示



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

视频

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

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

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

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

阅读文章 ➜

Delta德尔塔法

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

从这个属性中，我们可以得到（这是极值模型中使用的另一个参数化）的正态性，或者在任何四分位数上。我们运行一些模拟，再一次检查联合正态性。

> 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)

随时关注您喜欢的主题

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

   

``````
> n=5000

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

轮廓似然( Profile Likelihood )

另一个有趣的方法是Profile 似然函数的概念。因为尾部指数，在这里是辅助参数。
这可以用来推导出置信区间。在GPD的情况下，对于每个，我们必须找到一个最优的。我们计算Profile 似然函数，即。而我们可以计算出这个轮廓似然的最大值。

一般来说，这个两阶段的优化与（全局）最大似然是不等价的，计算结果如下

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

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

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

我们得到结果和最大似然估计的相似。我们可以用这种方法来计算置信区间，在图表上将其可视化

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

竖线是参数95%置信区间的下限和上限。

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。在此对他对本文所作的贡献表示诚挚感谢，他在上海财经大学完成了统计学专业的硕士学位，专注人工智能领域。擅长Python.Matlab仿真、视觉处理、神经网络、数据分析。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

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

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

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

随时关注您喜欢的主题

一般来说，这个两阶段的优化与（全局）最大似然是不等价的，计算结果如下

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

相关文章

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

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

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

随时关注您喜欢的主题

一般来说，这个两阶段的优化与（全局）最大似然是不等价的，计算结果如下

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

相关文章

关注有关新文章的微信公众号