R 中最大似然指数函数的参数

Question

我有一个样本数据，我正在尝试获取基于最大似然计算的双参数指数函数的参数。

我的样本：

sample = c(136.5,150,94.1,127.6,77.2,136.1,83.4,75.6,92.7,106.5,95.9,112.1,80.7,90.4,143.7,152.7,113.3,143.9,87.9,85.2,117.2,193,153.7,84.7,97.3,140.3,80,103.6,72.6,90.7,52.6,52.8)

我的主要目标是使用指数的 cdf 或 quantile 来获得最大似然，就像这样：

GEV 示例：

library(nsRFA)

parameters <- ML_estimation(sample, dist = "GEV")
p = c(0.1,0.066667,0.05,0.04,0.033333,0.02,0.01,0.005,0.002,0.001,0.0002,0.0001)
q = invF.GEV(1-p, parameters[1], parameters[2], parameters[3]); q
> 149.4 158.8 165.2 170 173.9 184.3 197.6 210 225.4 236.2 258.9 267.7

Answer 1

双参数指数函数是下端点在xi的指数函数。寻找具有如此尖锐边界点的分布的 MLE 有点特殊：边界的 MLE 等于数据集 中观察到的最小值（参见例如 this CrossValidated question).这使得双参数指数的 MLE 等同于 x-xmin.

的指数分布的 MLE

所以 xi 的 MLE 是

print(xi <- min(sample))

与MASS::fitdistr:

(m0 <- MASS::fitdistr(sample-xi, "exponential"))
      rate    
  0.018382353 
 (0.003249572)

与 bbmle:

m1 <- bbmle::mle2(y-xi~dexp(lambda), 
      data=data.frame(y=sample), start=list(lambda=1),
      method="Brent", lower=0.001, upper=100)
broom::tidy(m1, conf.int=TRUE, conf.method="profile")

  term   estimate std.error statistic      p.value conf.low conf.high
  <chr>     <dbl>     <dbl>     <dbl>        <dbl>    <dbl>     <dbl>
1 lambda   0.0184   0.00325      5.66 0.0000000154   0.0127    0.0255

解析解：

1/mean(sample-xi)
## [1] 0.01838235

如果您想要一个提供移位和缩放参数的简单函数（显然由您的替代软件提供）：

est_twoexp <- function(x) {
   xi <- min(x)
   c(xi = xi, scale = mean(sample-xi))
}
est_twoexp(sample)
##    xi scale 
##  52.6  54.4

est_twoexp <- function(x) {
   xi <- min(x)
   c(xi = xi, scale = mean(sample-xi))
}
est_twoexp(sample)
##    xi scale 
##  52.6  54.4

绘图：

library(nsRFA)
ee <- ecdf(sample)
inv_ecdf <- approxfun(seq(0, 1, length=length(sample)),
                      sort(sample),
                      method="constant")
## homemade quantile function
q2exp <- function(p, xi, scale) {
  xi + qexp(p, rate=1/scale)
}

curve(inv_ecdf(x), from=0, to=1, type="s", ylim=c(40,250))
with(as.list(est_twoexp(sample)),
     curve( invF.exp (x, xi, scale), col=2, lwd=2, add=TRUE))
with(as.list(est_twoexp(sample)),
     curve( q2exp(x, xi, scale), col=4, lwd=2, lty= 2, add=TRUE))

glm 和 family=Gamma 不起作用，因为它不允许零值（在 Gamma 分布的一般族中，x==0 只有正的有限密度指数分布）

R 中最大似然指数函数的参数

Parameters for Exponential function with maximum likelihood in R

r

exponential