R 中的 Weibull 参数估计，同时考虑 X（时间）和 Y（累积观察）

Posted 2023-02-19

【中文标题】R 中的 Weibull 参数估计，同时考虑 X（时间）和 Y（累积观察）

【英文标题】：Weibull parameter estimation in R, considering both X (time) and Y (cumulative observation)

【发布时间】：2021-06-06 06:36:09

【问题描述】：

早上好，

我正在尝试从以下数据集中计算 R 中的比例和形状 Weibull 参数

但是，X 数据是随时间不均匀抽取的样本，即它们不是均匀分布的。

X - 以天为单位的时间，可能从 0 到 14 不等 Y - 累积释放，从 0 到 100% 变化

我想请求您支持从该数据集中计算 Weibull 参数，因为根据我的研究，像 fitdist(dat, "weibull") 这样的函数似乎不适合 dat 是一个数字向量，我似乎无法将时间 (X) 对累积释放的影响包括在内。

X <- c(0.04,0.7,1,3,6,9)
Y <- c(12,46,48,67,87,93)

亲切的问候， 保罗

【问题讨论】：

fitdist 在这里不合适，因为它通过从 data 的随机样本中估计分布参数来工作。您可以从 CDF 中选择一些点。你需要做一些代数来求解相应的方程。 This 可能会有所帮助。

【参考方案1】：

您可以使用损失函数解决此问题，并注意 Weibull 分布的 CDF 由 F(x) = 1 - exp(-(x/lambda)^k) 给出。

library(tidyverse)

loss <- function(params) 
  x <- c(0.04,0.7,1,3,6,9)
  y <- c(12,46,48,67,87,93)/100
  print(paste0("lambda: ", params[1], "; k: ", params[2]))
  fitted <- 1 - exp(-(x/params[1])**params[2])
  loss <- sum((fitted - y)*(fitted - y))
  return(loss)


nlm(loss, c(1, 1))

tibble(
  x=rep(c(0.04,0.7,1,3,6,9)),
  y=c(12,46,48,67,87,93)/100,
  fitted=1 - exp(-(c(0.04,0.7,1,3,6,9)/1.95595756931587)**0.556026917419405)
) %>% 
  ggplot() +
  geom_point(aes(x=x, y=y)) +
  geom_line(aes(x=x, y=fitted))

给予

[1] "lambda: 1; k: 1"
[1] "lambda: 1; k: 1"
[1] "lambda: 1.000001; k: 1"
[1] "lambda: 1; k: 1.000001"
[1] "lambda: 1.22355098897031; k: 0.891679447194814"
[1] "lambda: 1.2235522125213; k: 0.891679447194814"
<lines deleted>
[1] "lambda: 1.95595952527344; k: 0.556025917419405"
[1] "lambda: 1.95595756931587; k: 0.556026917419405"

【讨论】：

亲爱的 Limey，非常感谢您的帮助，我能够实现它。亲切的问候，保罗

【参考方案2】：

模型可以适配nls。但y 的值必须在区间 [0, 1] 内。

df1 <- data.frame(x = X, y = Y/100)

fit <- nls(
  y ~ pweibull(x, shape = shape, scale = scale), 
  data = df1,
  start = list(shape = 0.1, scale = 0.1)
)
fit
#Nonlinear regression model
#  model: y ~ pweibull(x, shape = shape, scale = scale)
#   data: df1
#shape scale 
#0.556 1.956 
# residual sum-of-squares: 0.004963
#
#Number of iterations to convergence: 6 
#Achieved convergence tolerance: 6.215e-07

现在根据这个模型绘制一个图表。该模型似乎适合数据。

xnew <- seq(min(X), max(X), length.out = 100)
ynew <- predict(fit, newdata = data.frame(x = xnew))
newd <- data.frame(x = xnew, y = ynew)

plot(y ~ x, df1)
lines(y ~ x, newd, col = "red", lty = "dotted")

数据

X <- c(0.04,0.7,1,3,6,9)
Y <- c(12,46,48,67,87,93)

【讨论】：

快！你只赢了我几秒钟... ;)

@Limey 是的，刚刚注意到。估计与你的一致。

哈哈哈！是的，这总是一种解脱！我认为nls 比nlm 更整洁。我不知道如何从nlm 中获取参数估计值。 （第一次使用！）因此我的 hacky print...

@Limey 实际上，我喜欢打印的想法，带有参数verbose...

亲爱的Rui，非常感谢您的帮助，我能够实现它。亲切的问候，保罗