Fitting the the two-parameter Weibull Distribution

Introduction

This vignette presents the fitWD() function, which is used to estimate the parameters of the two-parameter Weibull Distribution (WD) using several estimation methods, including maximum likelihood (ML), least squares (LS), weighted least squares (WLS), and maximum product of spacings (MPS).

A random variable \(X\) is said to follow a two-parameter WD its cumulative distribution function (CDF) and probability density function (PDF) are given by \[ F(x) = 1- e^{-a x^b}, \] and \[ f(x) = a b x^{b-1} e^{-a x^b}, \]
where \(x>0\), \(a>0\) is the scale parameter and \(b > 0\) is a shape parameter.

Estimation methods

Let a random sample from \(X \sim WD(a, b)\) are observed as \(x_i, \; i=1,\dots, n\). For all estimation methods, optimization is performed using stats::optim.

The MLEs are obtained by maximizing the log-likelihood function: \[ \ell(a,b; \mathbf{x}) = \sum_{i=1}^{n} \log \big( a b x_i^{b-1} \big) -a \sum_{i=1}^{n} x_i^b. \]

The LSE is obtained by minimizing the squared differences between the theoretical and empirical CDF values at the ordered sample points.

Let \(x_{(1)}<x_{(2)}< \dots < x_{(n)}\) denote the ordered samples. The objective function is: \[ Q(a,b) = \sum_{i=1}^n \bigg( F(x_{(i)};a,b) - \frac{i-0.3}{n+0.4} \bigg)^2. \]

The WLSE extends LSE by introducing weights: \[ Q^W(a,b) = \sum_{i=1}^n w_i \bigg( F(x_{(i)};a,b) - \frac{i-0.3}{n+0.4} \bigg)^2, \] where \(w_i = \frac{(n+1)^2(n+2)}{i(n-i+1)}, \; i=1,\dots,n.\)

The MPS estimates are obtained by maximizing: \[ M(a,b) = \sum_{i=1}^{n+1} \frac{1}{n+1} \log \big( F(x_{(i)};a,b) - F(x_{(i-1)};a,b) \big), \] where \(F(x_{(0)};a,b) = 0\) and \(F(x_{(n+1)};a,b) = 1\).

Using fitWD()

To illustrate the use of fitWD(), we consider a simulated example.

# Load the package
library(SSReliabilityClaytonMWD)

# generate data from MWD(a, b, lambda) with lambda = 0 reduces two-parameter Weibull distribution
n <- 50
a <- 0.75; b <- 1.25; lambda <- 0
set.seed(123)
dat <- rMweibull(n, a, b, lambda)
# random initial points 
init <- runif(2)

MLE of the parameters:

# Fit WD to `dat`
fit.mle <- fitWD(data = dat, est.method = "mle", opt.method = "L-BFGS-B", starts = init,
                  lower = c(1e-05,1e-05), upper = c(Inf,Inf), hessian = FALSE )
fit.mle$estimates
#>        a        b 
#> 0.675279 1.283785

LSE of the parameters:

fit.lse <- fitWD(data = dat, est.method = "lse", opt.method = "L-BFGS-B", starts = init,
                   lower = c(1e-05,1e-05), upper = c(Inf,Inf), hessian = FALSE )
fit.lse$estimates
#>         a         b 
#> 0.7012337 1.2147202

WLSE of the parameters:

fit.wlse <- fitWD(data = dat, est.method = "wlse", opt.method = "L-BFGS-B", starts = init,
                    lower = c(1e-05,1e-05), upper = c(Inf,Inf), hessian = FALSE )
fit.wlse$estimates
#>        a        b 
#> 0.690969 1.236661

MPS of the parameters:

fit.mps <- fitWD(data = dat, est.method = "mps", opt.method = "L-BFGS-B", starts = init,
                   lower = c(1e-05,1e-05), upper = c(Inf,Inf), hessian = FALSE )
fit.mps$estimates
#>         a         b 
#> 0.6848478 1.1977846