极值理论（II 型和 III 型案例）优化问题（使用 R）

Question

我正在 R 中做一个关于极值理论的项目，我必须为 3 个分布中的每一个分布的参数编写我自己的最大似然估计。

极值分布的 mle 没有解析解，它们必须在数值上近似。 Gumbel 的情况很简单，因为我可以自己获取一个参数并使用结果在数值上优化另一个参数。

对于 Fréchet 和 Weibull 情况，您必须同时优化多个参数。我正在尝试在 R 中使用 optim() 函数，但得到的结果非常奇怪。我不确定这是我的代码有问题，我对 optim() 的理解有问题，还是理论本身有问题（或这 3 者的某种组合）。

无论如何，这是我的代码：

library(evd) #using this for frechet random variables

#generate maxima from 2 distributions: 1st should converge to Frechet, 2nd to Weibull
Mfrechet<-function (m) {
  d1<-matrix(rfrechet(1000000),nrow=m)
  apply(d1,1,max)/(1000000/m)
} 

Munif <- function (m) {
  d1<-matrix(runif(10000000),nrow=m)
  n<-ncol(d1)
    Xn<-apply(d1,1,max)
  n*Xn -n +1
}

x1<-Mfrechet(10000)
x2<-Munif(10000)

#negative log-likelihood for frechet
lf.fr <- function(y) {
  x<-x1
  n<-length(x1)
  a<-y[1]
  b<-y[2]
  g<-y[3]
  -(n*log(g)+n*g*log(b)-(g+1)*sum(log(x-a),na.rm=T)-sum((b/log(x-a))^g,na.rm=T))
}

#negative log likelihood function for weibull
lf.weibull <- function(y) {
  x<-x2
  n<-length(x2)
  a<-y[1]
  b<-y[2]
  g<-y[3]
  -(n*log(g)+n*g*log(b)+(g-1)*sum(log(a-x),na.rm=T)-sum(((a-x)/b)^g,na.rm=T))
}

optim(c(1,1,1),lf.fr)
optim(c(1,1,1),lf.weibull)

结果相去甚远，尤其是与 extRemes 包中的 fevd 相比：

library(extRemes)
fevd(x1)
fevd(x2)

我尝试更改 optim 的起始值，但它给了我截然不同的结果，但并不接近正确的结果。我用 getAnywhere 查看了不同包中的大多数 gev 拟合，它们大多也使用 optim，但它们很难分类。

谁能告诉我哪里出错了？希望有一种不太复杂的修复方法...

Answer 1

您的问题是two-fold。首先，GEV 分布具有取决于参数的支持。所以如果你知道参数，计算支持度是一件简单的事情，但如果你在优化器中搜索参数，很容易找到与你的数据不一致的组合。您的 loglik 函数需要处理该问题。

其次，GEV对应Weibull还是Frechet取决于形状参数xi（xi > 0表示Frechet；xi < 0表示Weibull）。所以你需要使用一个优化器，它允许你设置参数的下限和上限。如果您使用 optim(...)，您唯一的选择是“L-BFGS-B”。不幸的是，这种方法有些脆弱（参见 this post）。在同一个 post 中确定的替代方案是 dfoptim 包中的 nmkb(...)。

这是一个工作示例，基于 Wikipedia 中定义的 GEV 分布参数化。

library(evd)
library(dfoptim)
#
ll = \(par, x) {
  mu <- par[1]; sigma <- par[2]; xi <- par[3]
  s <- (x - mu)/sigma
  t <- (1 + xi*s)^(-1/xi)
  f <- t^(xi+1) * exp(-t) / sigma
  if(any(is.na(f)) | any(f <= 0)) return(1e6)
  return (-sum(log(f)))
}
par    <- c(mu=0, sigma=.5, xi=-0.1)
result <- nmkb(par, ll, lower = c(-Inf, 0, -Inf), upper = c(Inf, Inf, 0), x=x2)
result$par
## [1] -0.001692221  0.998771474 -0.997121661
fevd(x2)$results$par
##     location        scale        shape 
## -0.001711532  0.998739426 -0.997070437
#
par    <- c(mu=0, sigma=0.5, xi=+0.1)
result <- nmkb(par, ll, lower = c(-Inf, 0, 0), upper = c(Inf, Inf, Inf), x=x1)
result$par
## [1] 0.9997442 1.0003228 1.0004027
fevd(x1)$results$r
##  location     scale     shape 
## 0.9999759 1.0005464 1.0003256

请注意，您没有在 post 中设置随机数生成器的种子，因此您的 x1 和 x2 会有所不同。我用了 set.seed(1).

极值理论（II 型和 III 型案例）优化问题（使用 R）

Extreme Value Theory (Type II&III case) problem with optimization (using R)

statistics

optimization

r

log-likelihood