Reputation: 11
Good morning,
I am trying to calculate the scale and shape Weibull parameters in R, from the following dataset
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
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"
Upvotes: 1
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")
X <- c(0.04,0.7,1,3,6,9)
Y <- c(12,46,48,67,87,93)
Upvotes: 1