威布尔分布三参数

Question

我正在尝试使用最大似然法获取 Weibull 函数的三个参数。

我正在使用树木直径来创建林中树木的直径分布。但是，当我尝试获取它们时，出现错误：

警告：产生 NaN
警告：产生 NaN
警告：optim(c(gamma = 1, beta = 5, alpha = 10), ll.w3, data = DBH, 中产生NaNsError，： 非有限有限差分值 [1]

但是，这些错误仅在我使用我的数据集时发生；当我使用随机值时，代码完美运行，没有任何错误。

我正在使用书中的代码：“林业和环境数据的生物测定”。

使用随机值的代码：

set.seed(1)
diameters <- rnorm(30, mean = 50, sd = 6)
diameters


dweibull3 <- function(x, gamma, beta, alpha) {
 (gamma/beta)*((x - alpha)/beta)^(gamma - 1) *
  (exp(-((x - alpha)/beta)^gamma)) }

ll.w3 <- function(p, data)
   sum(log(dweibull3(data, p[1], p[2], p[3])))

mle.w3.nm <- optim(c(gamma = 1, beta = 5, alpha = 10),
                   ll.w3,
                   data =diameters,
                   hessian = TRUE,
                   control = list(fnscale = -1))
mle.w3.nm$par

结果：

  gamma      beta     alpha 
 6.843394 32.530438 20.122516 

---

使用我的数据集进行编码：

DBH <- my_datafile$DBH

mean_DBH <-mean(DBH)
sd_DBH <-sd(DBH)
DBH <- data.frame(DBH=DBH,Mean=mean_DBH, SD=sd_DBH)


dweibull3 <- function(x, gamma, beta, alpha) {
 (gamma/beta)*((x - alpha)/beta)^(gamma - 1) *
  (exp(-((x - alpha)/beta)^gamma)) }

ll.w3 <- function(p, data)
    sum(log(dweibull3(data, p[1], p[2], p[3])))

mle.w3.nm <- optim(c(gamma = 1, beta = 5, alpha = 10),
                   ll.w3,
                   data =DBH,
                   hessian = TRUE,
                   control = list(fnscale = -1))
mle.w3.nm$par

警告：产生 NaN
警告：产生了 NaN警告：产生了 NaN
optim(c(gamma = 1,beta = 5, alpha = 10), ll.w3, data = DBH, 中的错误：非有限有限差分值 [1]

Answer 1

威布尔尺度参数（此处为

beta

）的最大似然估计可以用形状参数（此处为gamma）以封闭形式来

表达

。这使我们能够仅使用两个参数来执行优化。下面的函数

w3mle

使用此策略来查找 3 参数 Weibull 的 MLE。

w3mle <- function(data) {
  n <- length(data)
  # initialize alpha to be near min(data)
  # initialize gamma parameter using quantile matching
  data <- sort(data)
  data0 <- data - data[1]
  q <- unname(quantile(data0 - data[1] + data[2], c(0.25, 0.75)))
  # initial parameters (search in log-space)
  p <- c(log(log1p(-0.25)/log1p(-0.75))/log(q[1]/q[2]), log(data0[2]))
  # objective function
  f <- function(p) {
    p <- exp(p)
    x <- data0 + p[2]
    beta <- (sum(x^p[1])/n)^(1/p[1])
    -sum(dweibull(x, p[1], beta, TRUE))
  }
  # optimization
  p <- exp(optim(p, f)$par)
  p[2] <- data[1] - p[2]
  c(gamma = p[1], beta = (sum((data - p[2])^p[1])/n)^(1/p[1]), alpha = p[2])
}

测试：

set.seed(1)
diameters <- rnorm(30, mean = 50, sd = 6)
(p <- w3mle(diameters))
#>     gamma      beta     alpha 
#>  19.15497  86.44453 -33.55127

根据问题中找到的参数集检查此解决方案：

p2 <- c(gamma = 6.843394, beta = 32.530438, alpha = 20.122516)
sum(dweibull(diameters - p[3], p[1], p[2], TRUE))
#> [1] -91.72539
sum(dweibull(diameters - p2[3], p2[1], p2[2], TRUE))
#> [1] -91.99055

p

是比

p2

稍好的估计。

现在尝试一下

DBH

：

w3mle(DBH)
#>    gamma     beta    alpha 
#> 3.960256 3.504212 6.493401