The Power Generalised Weibull distribution

The PGW pdf, survival function, and hazard functions of the PGW are given below [see also PGW]:

\[ f(t \mid \eta,\nu,\delta) = \dfrac{\nu}{\delta \eta^\nu}t^{\nu-1} \left[ 1 + \left(\dfrac{t}{\eta}\right)^\nu\right]^{\left(\frac{1}{\delta}-1\right)} \exp\left\{ 1- \left[ 1 + \left(\dfrac{t}{\eta}\right)^\nu\right]^{\frac{1}{\delta}} \right\}, \]

\[ S(t \mid \eta,\nu,\delta) = \exp\left\{ 1- \left[ 1 + \left(\dfrac{t}{\eta}\right)^\nu\right]^{\frac{1}{\delta}} \right\}, \]

\[ h(t \mid \eta,\nu,\delta) = \dfrac{\nu}{\delta \eta^\nu}t^{\nu-1} \left[ 1 + \left(\dfrac{t}{\eta}\right)^\nu\right]^{\left(\frac{1}{\delta}-1\right)}. \] where \(\eta>0\) is a scale parameter, and \(\nu,\delta >0\) are shape parameters.

Total variation distance between PGW and Weibull distributions

Alvares and Rubio (2021) proved the following result.

Theorem. The Total Variation distance between the PGW distribution with parameters \((\eta,\nu,\delta)\) and a Weibull distribution with parameters \((\eta,\nu)\) only depends on the parameter \(\delta\). This is:

\[ d_{TV}(f_{PGW},f_{W} \mid \eta, \nu, \delta) = d_{TV}(f_{PGW},f_{W} \mid 1, 1, \delta). \]

R code

The following R code shows the TV distance between the PGW and Weibull distributions, as well as the corresponding BTV(1,1) prior. For more information about BTV priors, see Dette, Ley, and Rubio (2018). The best approximation is determined using least squares, leading to a Gamma prior with scale parameter \(1.83\) and shape parameter \(0.65\). The \(BTV(1,1)\) prior and its approximation are shown in the following Figure.

rm(list=ls())

# Probability Density Function PGW
dpgw <- function(t, sigma, nu, gamma, log = FALSE){
  val <- log(nu) - log(gamma) - nu*log(sigma) + (nu-1)*log(t) + 
    (1/gamma - 1)*log( 1 + (t/sigma)^nu ) + 
    ( 1 - ( 1 + (t/sigma)^nu )^(1/gamma) )
  out <- ifelse(log, val, exp(val))
  return(out)
}

# Total variation distance between a PGW and Weibull distributions as a function of the parameter gamma
nu = 1
dtv <- Vectorize(function(gamma){
  tempf <- Vectorize(function(t) 0.5*abs( dpgw(t,1,nu,gamma) -  dweibull(t, shape = nu, scale = 1)) )
  val <- integrate(tempf,0,Inf, subdivisions = 1000)$value
  return(val)
} ) 

curve(dtv,0.01,25, xlab =expression(delta), ylab = "TV Distance", lwd = 2, cex.axis = 1.5, cex.lab = 1.5, n = 500)

# Perturbation measure between a PGW and Weibull distributions as a function of the paramter gamma
nu = 1
mtv <- Vectorize(function(gamma){
  tempf <- Vectorize(function(t) 0.5*abs( dpgw(t,1,nu,gamma) -  dweibull(t, shape = nu, scale = 1)) )
  val <- integrate(tempf,0,Inf)$value
  val <- ifelse(gamma < 1,-val,val)
} ) 

curve(mtv,0.01,25, xlab =expression(delta), ylab = "Perturbation measure", lwd = 2, cex.axis = 1.5, cex.lab = 1.5, n = 500)

# BTV(1,1) prior for gamma
# https://doi.org/10.1111/sjos.12306
library(numDeriv)
prior <- Vectorize(function(gamma){
  val <- abs( grad(func = mtv, x = gamma ))
  return(val/1.869813)
})

curve(prior,0.05,5, n = 50)

# Best Gamma approximation using least squares
n = 50
vec <- seq(0.25,4,length.out = n)
p.vec <- prior(vec)
# Sum of Squared Errors
sse <- function(par){
  par0 <- exp(par)
  sse <- sum( ( p.vec - dgamma(vec, shape = par0[1],  scale = par0[2]) )^2 ) 
  return(sse)
}

OPT = optim(c(0,0), sse, control = list(maxit = 1000))

best.w <- Vectorize(function(t) dgamma(t,  shape = exp(OPT$par[1]), scale = exp(OPT$par[2])))

plot(vec, p.vec, type = "l", col = "blue", lwd = 2)
curve(best.w, 0.001,7, col = "red", add=T, n = 1000)

# Parameters
exp(OPT$par)
## [1] 0.645247 1.833319
best.w <- Vectorize(function(t) dgamma(t,  shape = 0.65, scale = 1.83))

n=100
vec <- seq(0.01,7,length.out = n)
vec <- vec[-15]
p.vec <- prior(vec)

plot(vec, p.vec, type = "l", col = "black", lwd = 2, xlab = ~delta, ylab = "Density",
     cex.lab = 1.5, cex.axis = 1.5)
curve(best.w, 0.001,7, col = "black", lwd = 2, lty = 2, add=T, n = 1000)
legend(4,0.9, c("BTV(1,1)","Approximation"),
       text.col = c("black","black"), col = c("black","black"), 
       lty = c(1,2),lwd=2, merge = TRUE, bg = "gray90")

References

Alvares, D., and F. J. Rubio. 2021. “A Tractable Bayesian Joint Model for Longitudinal and Survival Data.” Statistics in Medicine.

Dette, H., C. Ley, and F. J. Rubio. 2018. “Natural (Non-) Informative Priors for Skew-Symmetric Distributions.” Scandinavian Journal of Statistics 45 (2): 405–20. https://doi.org/10.1111/sjos.12306.