Paulo Francisco
Paulo Francisco

Reputation: 11

Weibull parameter estimation in R, considering both X (time) and Y (cumulative observation)

Good morning,

I am trying to calculate the scale and shape Weibull parameters in R, from the following dataset

enter image description here

However, the X data are samples pulled unevenly across time, i.e., they are not equally distributed.

I would like to ask for your support in calculating the Weibull parameters from this data set, as from what I was able to research, functions such as fitdist(dat, "weibull") do not seem to be approriatte as dat is a vector of numbers, and I seem unable to include the influence of time (X) over the cummulative release.

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

Kind regards, Paulo

Upvotes: 1

Views: 905

Answers (2)

Limey
Limey

Reputation: 12586

You can solve this problem using a loss function and noting that the CDF of the Weibull distributiuon is given by 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))

Giving

[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"

enter image description here

Upvotes: 1

Rui Barradas
Rui Barradas

Reputation: 76641

The model can be fitted with nls. But the y values must be in the interval [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

Now plot a graph based on this model. The model seems to fit the data.

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

enter image description here

Data

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

Upvotes: 1

Related Questions