# Modified Bivariate Weibull Distribution Allowing Instantaneous and Early
Failures

Sumangal Bhattacharya  Ishapathik Das<EMAIL_ADDRESS>Muralidharan
Kunnummal  Department of Mathematics, Indian Institute of Technology Tirupati,
Tirupati, India.  Department of Statistics, The M. S. University of Baroda,
Vadodara, India.

###### Abstract

In reliability and life data analysis, the Weibull distribution is widely used
to accommodate more data characteristics by changing the values of the
parameters. We frequently observe many zeros or close to zero data points in
reliability and life testing experiments. We call this phenomenon a nearly
instantaneous failure. Many researchers modified the commonly used univariate
parametric models such as exponential, gamma, Weibull, and log-normal
distributions to appropriately fit such data having instantaneous failure
observations. Researchers also find bivariate correlated life testing data
having many observations near a point $(x_{0},y_{0})$ while the remaining
observations follow some continuous distribution. This situation defines as
responses having early failures for such bivariate responses. In particular,
if the point is the origin [i.e., $(x_{0},y_{0})=(0,0)$], then we call the
situation a nearly instantaneous failure for the responses. Here, we propose a
modified bivariate Weibull distribution that allows early failure by combining
bivariate uniform distribution and bivariate Weibull distribution. The
bivariate Weibull distribution is constructed using a 2-dimensional copula,
assuming the marginal distributions as two parametric Weibull distributions.
We derive some properties of that modified bivariate Weibull distribution,
mainly the joint probability density function, the survival (reliability)
function, and the hazard (failure rate) function. The model’s unknown
parameters are estimated using the Maximum Likelihood Estimation (MLE)
technique combined with a machine learning clustering algorithm. Numerical
examples are provided using simulated data to illustrate and test the
performance of the proposed methodologies. The method is also applied to real
data and compared with existing approaches to model such data in the
literature.

###### keywords:

Bivariate models, copula, DBSCAN clustering algorithm, Instantaneous failure,
Hazard function, Survival function, Weibull distribution.

##  1 Introduction

In reliability analysis and life data analysis, Weibull distribution is used
mostly because of its versatility. In literature, extension and generalization
of the Weibull distribution are given in Pham and Lai, (2007), Murthy et al.,
(2004), El-Gohary and El-Morshedy, (2015). In life testing experiments, mainly
in electronic parts and clinical trials, many items fail instantaneously or
fail early due to faulty construction or inferior quality or due to no
response of the treatments. We record the observed lifetime is close to zero.
To model life data containing such nearly instantaneous failure, traditional
parametric distributions, for example, Weibull distribution, Pareto
distribution, exponential distribution, gamma distribution, log-normal
distribution, etc., are not accurate in practice as such failures in
experiment usually discard the assumption of a uni-modal distribution. These
situations can be modeled by modifying the traditionally used parametric
distributions. Muralidharan and Lathika, (2006), modified the two parameters
Weibull distribution by mixing two-parameter Weibull distribution and a
singular distribution at zero. Lai et al., (2007) modified the two parameters
Weibull distribution by considering the mixture of the uniform distribution
and two parameters Weibull distribution. Muralidharan and Khabia, (2011),
model this instantaneous or early failures data by considering the mixture of
the uniform distribution and the Pareto distribution. To estimate the model
parameters, maximum likelihood estimation is used mainly. Muralidharan et al.,
(2011), estimated the model parameters using Bayesian approach. Cheng and
Sheu, (2016) proposed M-estimation method and derived the closed form of the
shape parameter estimation. Using simulation study they showed that,
M-estimation performs better than the MLE for samples with early failure. To
model data with inliers (instantaneous or early failure), non-standard mixture
of distributions, with degeneracy occurring at zero and one, and a probability
distribution for positive observation is used mainly. Muralidharan and
Bavagosai, (2018) estimates the parameters of such model on the basis of type-
II censored sample from a Weibull distribution with discrete mass at zero and
one. To model instantaneous failure data from clinical trials and
longitudinal, Withana Gamage et al., (2020) proposed two methods to estimate
the parameters and used expectation-maximization(EM) algorithm. Muralidharan
and Bavagosai, (2021) provided a method for testing hypothesis to find
statistical significance of the single and multiple inliers in samples.

Many data in real life are bivariate in nature and correlated, which contains
huge observation near a particular point $(x_{0},y_{0})$ and other responses
follow some continuous distribution; we define this as early failure data. If
that particular point is origin, i.e. $(x_{0},y_{0})=(0,0)$, we call it as
nearly instantaneous failure data. Bivariate Weibull distribution using FGM
copula and there some properties are proposed in literature Almetwally et al.,
(2020), Peres et al., (2018). However, we did not find any modeling technique
to fit such correlated bivariate data with instantaneous or nearly
instantaneous failure observations adequately in the literature. Here, we
propose a modified bivariate Weibull distribution allowing early failure by
considering the mixture of the bivariate uniform distribution and bivariate
Weibull distribution. The bivariate Weibull distribution is obtained using a
copula and two parameters Weibull distribution as marginal distribution. The
forms of the survival function and hazard function are obtained by using a
generalization of the Farlie-Gumbel-Morgenstern (GFGM) copula given by Lai and
Xie, (2000). Furthermore, we consider the model to deal with the presence of
nearly instantaneous failure $[(x_{0},y_{0})=(0,0)]$ in the responses and
obtained the form of the hazard function and survival function using the GFGM
copula. To illustrate the proposed model, in the simulation study, we consider
the Gaussian copula [Nelsen, (2007)] to construct the modified bivariate
Weibull distribution. To estimate the parameters, we use the method of maximum
likelihood estimation (MLE) and density-based spatial clustering of
applications with noise (DBSCAN) clustering given by Ester et al., (1996).

The remainder of the article is organized as follows. We provide a brief
review of copula functions and Sklar’s theorem for constructing multivariate
distribution functions using the known marginals in Section 2. In Section 3,
we construct a bivariate Weibull distribution function using a copula having
two parameters Weibull distributions as marginals. We propose a modified
Weibull distribution for early and instantaneous failure bivariate
observations along with their hazard functions and the survival functions in
Sections 4 and 5. Estimation of unknown parameters using the MLE combined with
DBSCAN clustering process is described in Section 6. In Section 7, we provide
an algorithm to simulate data from modified bivariate Weibull distribution.
Numerical examples with simulated data are given in Section 8.1, and the
proposed method is applied to real data in 8.2. The concluding remarks are
given in Section 9.

##  2 Copulas and Sklar’s Theorem

In this section, copula function and how it can be used to form multivariate
joint probability mass function from its marginals are discussed. For more
information about copula we refer Nelsen, (2007), Joe, (2014) and Song and
Song, (2007).

######  Definition 1.

A d-dimensional copula C is a function from $[0,1]^{d}$ to $[0,1]$ satisfying
the following conditions:

  1. 1.

$C(1,\ldots,1,y_{j},1,\ldots,1)=y_{j},\forall\ j=1,2,\ldots,d$ where
$(1,\ldots,1,y_{j},1,\ldots,1)\in[0,1]^{d}$,

  2. 2.

$C(y_{1},y_{2},\ldots,y_{d})=0$ if at least one $y_{j}=0$ for
$j=1,2,\ldots,d$,

  3. 3.

For any $\textbf{y}_{1}=(y^{(1)}_{1},y^{(1)}_{2},\cdots,y^{(1)}_{d})$,
$\textbf{y}_{2}=(y^{(2)}_{1},y^{(2)}_{2},\cdots,y^{(2)}_{d})$ belong to
$[0,1]^{d}$ with $y^{(1)}_{j}\leq y^{(2)}_{j}$, for all $j=1,2,\ldots,d$:

$\displaystyle\sum_{i_{1}=1}^{2}\sum_{i_{2}=1}^{2}\ldots\sum_{i_{d}=1}^{2}(-1)^{i_{1}+i_{2}+\ldots+i_{d}}C(y^{(i_{1})}_{1},y^{(i_{2})}_{2},\ldots,y^{(i_{d})}_{d})\geq
0.$ (2.1)

In copula theory, Sklar’s theorem [Sklar, (1959)] is a fundamental result and
we can use it to form a joint distribution from its marginal distributions.

######  Theorem 1.

Let $F_{1},F_{2},\ldots,F_{d}$ are the marginal distribution functions for the
random variables $Y_{1},Y_{2},\ldots,Y_{d}$ and their joint cumulative
distribution function be F, then the following are true:

  1. 1.

$\exists$ a d-dimensional copula C s.t., $\forall$
$(y_{1},y_{2},\ldots,y_{d})\in\mathbb{R}^{d}$,

$\displaystyle
F(y_{1},y_{2},\ldots,y_{d})=C(F_{1}(y_{1}),F_{2}(y_{2}),\ldots,F_{d}(y_{d})),$
(2.2)

  2. 2.

If $Y_{1},Y_{2},\ldots,Y_{d}$ are continuous then the copula C is unique.
Otherwise, C can be uniquely determined on a d-dimensional rectangle
$Range(F_{1})\times Range(F_{2})\times\ldots\times Range(F_{d}).$

A copula is a multivariate probability distribution with uniform marginal
distributions on the interval $[0,1]$. Copulas are popular because of their
ability to model dependence among $Y_{1},Y_{2}\cdots,Y_{d}$. Here is some
examples of copula

Example 1.  The multivariate Gaussian copula is given by the function:

$\displaystyle
C(y_{1},y_{2},\ldots,y_{d}|R)=\bm{\Phi}_{R}(\Phi^{-1}(y_{1}),\Phi^{-1}(y_{2}),\ldots,\Phi^{-1}(y_{d})),$
(2.3)

where $(y_{1},y_{2},\ldots,y_{d})\in[0,1]^{d}$, $\Phi^{-1}$ is the inverse of
cumulative distribution function of the standard Gaussian distribution and
$\bm{\Phi}_{R}$ is the joint CDF of a standard multivariate Gaussian
distribution with correlation matrix $R$.

The Gaussian copula density (Arbenz, (2013)) is given by

$\displaystyle c(y_{1},y_{2},\ldots,y_{d}\mid R)=\frac{1}{\sqrt{\mid
R\mid}}exp\left\\{-\frac{1}{2}\textbf{U}^{T}\times\left(R^{-1}-I_{d}\right)\times\textbf{U}\right\\},$
(2.4)

where
$\textbf{U}=\left[\Phi^{-1}(y_{1}),\Phi^{-1}(y_{2}),\ldots,\Phi^{-1}(y_{d})\right]^{T}$,
$(y_{1},y_{2},\ldots,y_{d})\in[0,1]^{d}$, $\Phi$ is the univariate standard
normal distribution and $R$ is the correlation matrix of the standard
multivariate Gaussian distribution and $I_{d}$ be the $d\times d$ identity
matrix.

Example 2.  A generalization of the bivariate Farlie-Gumbel-Morgenstern (GFGM)
copula Lai and Xie, (2000) is given by

$C(x,y)=xy+\rho x^{b}y^{b}(1-x)^{a}(1-y)^{a}$ (2.5)

where $(x,y)\in[0,1]^{2}$, $a\geq 1,\ b\geq 1$ and $-1\leq\rho\leq 1$.

Therefore, the density function of the GFGM copula is given by

$\begin{split}c(x,y)=&1+b^{2}\rho x^{b-1}y^{b-1}(1-x)^{a}(1-y)^{a}-ab\rho
x^{b-1}y^{b}(1-x)^{a}(1-y)^{a-1}-\\\ &ab\rho
x^{b}y^{b-1}(1-x)^{a-1}(1-y)^{a}+a^{2}\rho
x^{b}y^{b}(1-x)^{a-1}(1-y)^{a-1}\end{split}$ (2.6)

where $(x,y)\in[0,1]^{2}$, $a\geq 1,\ b\geq 1$ and $-1\leq\rho\leq 1$.

##  3 Bivariate Weibull Distribution

Let us recall the definition of two parameters Wibull distribution. We denote
$W\sim Wb(\alpha,\beta)$, if a random variable $W$ follows Weibull
distribtuion with parameters $(\alpha,\beta)$ with the CDF given by

$F_{W}(x)=1-\exp\left[-\left(\dfrac{x}{\beta}\right)^{\alpha}\right],\quad\quad
x>0,\alpha>0,\beta>0$ (3.1)

and the PDF given by

$f_{W}(x)=\dfrac{\alpha}{\beta}\left(\dfrac{x}{\beta}\right)^{\alpha-1}\exp\left[-\left(\dfrac{x}{\beta}\right)^{\alpha}\right],\quad\quad
x>0,\ \alpha>0,\ \beta>0$ (3.2)

Let $X\sim Wb(\alpha_{1},\beta_{1})$ and $Y\sim Wb(\alpha_{2},\beta_{2})$
follow two parameter Weibull distribution with parameters
$(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ respectively with
$\alpha_{1},\beta_{1},\alpha_{2},\beta_{2}>0$. Therefore the joint CDF of the
bivariate Weibull distribution using a 2-dimensional copula $C$ and Sklar’s
theorem is given by

$F_{XY}(x,y)=C(F_{X}(x),F_{Y}(y))$ (3.3)

where $x>0$, $y>0$ and $F_{X},F_{Y}$ are the CDF of the random variables X and
Y respectively.

The joint probability density function (PDF) of the bivariate Weibull
distribution is given by

$f_{XY}(x,y)=f_{X}(x)f_{Y}(y)\dfrac{\partial^{2}C(F_{X}(x),F_{Y}(y))}{\partial
F_{X}\partial F_{Y}}$ (3.4)

where $x>0$, $y>0$ and $f_{X}(x),\ f_{Y}(y)$ are the PDF of the random
variables X and Y respectively.

Now, we consider the GFGM copula (2.6) to find the form of the joint CDF and
joint PDF of the bivariate Weibull distribution. For simplicity of the
expressions, we define notation for some commonly used term as

$\begin{split}A=&\left(\dfrac{x}{\beta_{1}}\right)^{\alpha_{1}}\\\
B=&\left(\dfrac{y}{\beta_{2}}\right)^{\alpha_{2}}\\\ C=&A+B\\\
D=&\exp\left(-(a-1)C\right)\left((a+b)\exp\left(-A\right)-a\right)\left((a+b)\exp\left(-B\right)-a\right).\end{split}$
(3.5)

Therefore, the joint CDF of the bivariate Weibull distribution using the GFGM
copula (2.6) is given by

$F_{XY}(x,y)=F_{X}(x)F_{Y}(y)+\rho
F_{X}(x)^{b}F_{Y}(y)^{b}(1-F_{X}(x))^{a}(1-F_{Y}(y))^{a}$ (3.6)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $a\geq 1$, $b\geq 1$, $-1\leq\rho\leq 1$ and $F_{X}(x)$,
$F_{Y}(y)$ are the CDF of $X$ and $Y$ respectively.
Therefore, the joint PDF $f_{XY}$ of the bivariate Weibull distribution is
given by

$\begin{split}f_{XY}(x,y)=&\dfrac{\alpha_{1}\alpha_{2}}{\beta_{1}\beta_{2}}A^{1-1/\alpha_{1}}B^{1-1/\alpha_{2}}\exp(-C)[1+\rho\left(1-\exp\left(-A\right)\right)^{(b-1)}\left(1-\exp\left(-B\right)\right)^{(b-1)}D]\end{split}$
(3.7)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $a\geq 1$, $b\geq 1$, $-1\leq\rho\leq 1$.

By taking restrictions on $a$ and $b$, we can get a simplified form of the
bivariate Weibull distribution functions.

  * 1.

If we choose $a=1,b=1$, the form of the joint CDF is

$\begin{split}F_{XY}(x,y)=&F_{X}(x)F_{Y}(y)[1+\rho(1-F_{X}(x))(1-F_{Y}(y))]\\\
=&\left[1-\exp\left(-A\right)\right]\left[1-\exp\left(-B\right)\right]\left[1+\rho\exp\left[-\left(C\right)\right]\right]\end{split}$
(3.8)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $0\leq\rho\leq 1$.
Therefore the corresponding joint PDF of the bivariate Weibull distribution is

$\begin{split}f_{XY}(x,y)=&\dfrac{\alpha_{1}\alpha_{2}}{\beta_{1}\beta_{2}}A^{1-1/\alpha_{1}}B^{1-1/\alpha_{2}}\exp(-C)\left[1+\rho\left(2\exp\left(-A\right)-1\right)\left(2\exp\left(-B\right)-1\right)\right]\end{split}$
(3.9)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $-1\leq\rho\leq 1$.

##  4 Modified Bivariate Weibull Distribution

Here we find the distribution function of the modified Weibull distribution
and the form of the survival function and hazard function.

###  4.1 PDF and CDF of modified bivariate Weibull distribution

We consider the modified bivariate model allowing the early failure, as a
mixture of the bivariate uniform distribution and the bivariate Weibull
distribution. Thus the joint PDF of the modified bivariate Weibull (MBW)
distribution is of the form

$\begin{split}f(x,y)&=pf_{1}(x,y)+qf_{2}(x,y)\\\ \end{split}$ (4.1)

where $p+q=1$, $0<p<1$, $f_{1}(x,y)$ be the bivariate uniform distribution on
the rectangle $[t_{0},t_{0}+d]\times[t_{0},t_{0}+d]$, given by

$f_{1}(x,y)=\left\\{\begin{array}[]{ll}\dfrac{1}{d^{2}},&\quad x_{0}\leq x\leq
x_{0}+d,\ y_{0}\leq y\leq y_{0}+d\\\
0,&\quad\text{Otherwise}\end{array}\right.$ (4.2)

for sufficiently small $d(>0)$, and $f_{2}(x,y)=f_{XY}(x,y)$ is the bivariate
Weibull distribution, given in (3.4).
Therefore, the joint CDF of the MBW distribution is,

$F(x,y)=pF_{1}(x,y)+qF_{2}(x,y)$ (4.3)

where $F_{1}$ and $F_{2}$ are joint CDF of the bivariate uniform distribution
and the bivariate Weibull distribution respectively.

###  4.2 Survival function and Hazard function

The survival function of the MBW distribution is given by

$\begin{split}R(x,y)&=F((X,Y)>(x,y))\\\
&=pF_{1}((X,Y)>(x,y))+qF_{2}((X,Y)>(x,y))\\\
&=pR_{1}(x,y)+qR_{2}(x,y)\end{split}$ (4.4)

where $R_{1}$ and $R_{2}$ are the survival functions of the bivariate uniform
distribution and the bivariate Weibull distribution respectively.
The hazard function of the MBW distribution is given by

$h(x,y)=\dfrac{\dfrac{\partial^{2}R(x,y)}{\partial x\partial
y}}{R(x,y)}=\dfrac{f(x,y)}{R(x,y)}.$ (4.5)

Now, the survival functions $R_{1}$, $R_{2}$ and the hazard functions $h_{1}$,
$h_{2}$ of the respective component distributions are given respectively

$R_{1}(x,y)=\left\\{\begin{array}[]{ll}1,&\quad x\leq x_{0},\ y\leq y_{0}\\\
\dfrac{x_{0}+d-x}{d},&\quad x_{0}\leq x\leq x_{0}+d,\ y\leq y_{0}\\\
\dfrac{y_{0}+d-y}{d},&\quad x\leq x_{0},\ y_{0}\leq y\leq y_{0}+d\\\
\dfrac{(x_{0}+d-x)(y_{0}+d-y)}{d^{2}},&\quad x_{0}\leq x\leq x_{0}+d,\
y_{0}\leq y\leq y_{0}+d\\\ 0,&\quad x_{0}+d\leq x\ \text{or}\ y_{0}+d\leq
y\end{array}\right.$ (4.6)

$h_{1}(x,y)=\left\\{\begin{array}[]{ll}\dfrac{1}{(x_{0}+d-x)(y_{0}+d-y)},&\quad
x_{0}\leq x<x_{0}+d,\ y_{0}\leq y<y_{0}+d\\\ \infty,&\quad x_{0}+d\leq x\
\text{or}\ y_{0}+d\leq y\\\ 0,&\quad\text{otherwise}\end{array}\right.$ (4.7)

$\begin{split}R_{2}(x,y)=&1-F_{X}(x)-F_{Y}(y)+F_{XY}(x,y)\\\
=&1-F_{X}(x)-F_{Y}(y)+C(F_{X}(x),F_{Y}(y))\end{split}$ (4.8)

Where $F_{X}$, $F_{Y}$, and $F_{XY}$ are the CDF of the random variables X, Y,
and the bivariate Weibull distribution, respectively, and $C$ is a bivariate
copula.

$h_{2}(x,y)=\dfrac{f_{X,Y}(x,y)}{R_{2}(x,y)}$ (4.9)

where the joint PDF of the bivariate Weibull distribution $f_{XY}(x,y)$ is
given in (3.4) and the survival function $R_{2}(x,y)$ is given by (4.8).

In particular, if we consider copula as GFGM copula, then the survival
function $R_{2}(x,y)$ and hazard function $h_{2}(x,y)$ of the bivariate
Weibull distribution is given by

$\begin{split}R_{2}(x,y)=&\exp\left(-C\right)\left[1+\rho\exp\left(-(a-1)C\right)\left(1-\exp\left(-A\right)\right)^{b}\left(1-\exp\left(-B\right)\right)^{b}\right]\end{split}$
(4.10)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $a\geq 1$, $b\geq 1$, $-1\leq\rho\leq 1$.

and

$\begin{split}h_{2}(x,y)=&\dfrac{\begin{split}\dfrac{\alpha_{1}\alpha_{2}}{\beta_{1}\beta_{2}}A^{1-1/\alpha_{1}}B^{1-1/\alpha_{2}}\exp(-C)\left[1+\rho\left(1-\exp\left(-A\right)\right)^{(b-1)}\left(1-\exp\left(-B\right)\right)^{(b-1)}D\right]\end{split}}{\begin{split}\exp\left(-C\right)\left[1+\rho\exp\left(-(a-1)C\right)\left(1-\exp\left(-A\right)\right)^{b}\left(1-\exp\left(-B\right)\right)^{b}\right]\end{split}}\end{split}$
(4.11)

where $x>0$, $y>0$, $\alpha_{1}>0$, $\beta_{1}>0$, $\alpha_{2}>0$,
$\beta_{2}>0$, $a\geq 1$, $b\geq 1$, $-1\leq\rho\leq 1$.

The hazard function for the MBW distribution (4.5) can be expressed as

$h(x,y)=\dfrac{f(x,y)}{R(x,y)}=w(x,y)h_{1}(x,y)+(1-w(x,y))h_{2}(x,y)$ (4.12)

where $w(x,y)=\dfrac{pR_{1}(x,y)}{R(x,y)}$ for all $x>0,\ y>0$.
Therefore,

$w(x,y)=\left\\{\begin{array}[]{ll}\dfrac{p}{R(x,y)},&\quad x\leq x_{0},\
y\leq y_{0}\\\ \dfrac{p(x_{0}+d-x)}{dR(x,y)},&\quad x_{0}\leq x\leq x_{0}+d,\
y\leq y_{0}\\\ \dfrac{p(y_{0}+d-y)}{dR(x,y)},&\quad x\leq x_{0},\ y_{0}\leq
y\leq y_{0}+d\\\ \dfrac{p(x_{0}+d-x)(y_{0}+d-y)}{d^{2}R(x,y)},&\quad x_{0}\leq
x\leq x_{0}+d,\ y_{0}\leq y\leq y_{0}+d\\\ 0,&\quad x_{0}+d\leq x\ \text{or}\
y_{0}+d\leq y\end{array}\right.$ (4.13)

Thus, the expression for $R(x,y)$, $f(x,y)$ and $h(x,y)$ are given by

$R(x,y)=\left\\{\begin{array}[]{ll}p+qR_{2}(x,y),&\quad x\leq x_{0},\ y\leq
y_{0}\\\ \dfrac{p(x_{0}+d-x)}{d}+qR_{2}(x,y),&\quad x_{0}\leq x\leq x_{0}+d,\
y\leq y_{0}\\\ \dfrac{p(y_{0}+d-y)}{d}+qR_{2}(x,y),&\quad x\leq x_{0},\
y_{0}\leq y\leq y_{0}+d\\\
\dfrac{p(x_{0}+d-x)(y_{0}+d-y)}{d^{2}}+qR_{2}(x,y),&\quad x_{0}\leq x\leq
x_{0}+d,\ y_{0}\leq y\leq y_{0}+d\\\ qR_{2}(x,y),&\quad x_{0}+d\leq x\
\text{or}\ y_{0}+d\leq y\end{array}\right.$ (4.14)
$f(x,y)=\left\\{\begin{array}[]{ll}\dfrac{p}{d^{2}}+qf_{XY}(x,y),&\quad
x_{0}\leq x\leq x_{0}+d,\ y_{0}\leq y\leq y_{0}+d\\\
qf_{XY}(x,y),&\quad\text{Otherwise}\end{array}\right.$ (4.15)

$h(x,y)=\left\\{\begin{array}[]{ll}\dfrac{qf_{XY}(x,y)}{p+qR_{2}(x,y)},&\quad
x\leq x_{0},\ y\leq y_{0}\\\
\dfrac{qdf_{XY}(x,y)}{p(x_{0}+d-x)+qdR_{2}(x,y)},&\quad x_{0}\leq x\leq
x_{0}+d,\ y\leq y_{0}\\\
\dfrac{qdf_{XY}(x,y)}{p(y_{0}+d-y)+qdR_{2}(x,y)},&\quad x\leq x_{0},\
y_{0}\leq y\leq y_{0}+d\\\
\dfrac{p+qd^{2}f_{XY}(x,y)}{p(x_{0}+d-x)(y_{0}+d-y)+qd^{2}R_{2}(x,y)},&\quad
x_{0}\leq x\leq x_{0}+d,\ y_{0}\leq y\leq y_{0}+d\\\
\dfrac{f_{XY}(x,y)}{R_{2}(x,y)},&\quad x_{0}+d\leq x\ \text{or}\ y_{0}+d\leq
y\end{array}\right.$ (4.16)

##  5 Nearly instantaneous failure case $\left[(x_{0},y_{0})=(0,0)\right]$

Consider $(x_{0},y_{0})=(0,0)$ as a special case of the modified model (4.1),
which is called as nearly instantaneous failure model. In this case the
survival functions and hazard function for the bivariate uniform distribution
are given by

$R_{1}(x,y)=\left\\{\begin{array}[]{ll}\dfrac{(d-x)(d-y)}{d^{2}},&\quad 0\leq
x\leq d,\ 0\leq y\leq d\\\ 0,&\quad d\leq x\ \text{or}\ d\leq
y\end{array}\right.$ (5.1)

and

$h_{1}(x,y)=\left\\{\begin{array}[]{ll}\dfrac{1}{(d-x)(d-y)},&\quad 0\leq
x\leq d,\ 0\leq y\leq d\\\ \infty,&\quad d\leq x\ \text{or}\ d\leq
y.\end{array}\right.$ (5.2)

The survival functions and hazard function for the modified bivariate Weibull
model with nearly instantaneous failure occurring uniformly over the rectangle
$[0,d]\times[0,d]$ is given by

$R(x,y)=\left\\{\begin{array}[]{ll}\dfrac{p(d-x)(d-y)}{d^{2}}+qR_{2}(x,y),&\quad
0\leq x\leq d,\ 0\leq y\leq d\\\ qR_{2}(x,y),&\quad d\leq x\ \text{or}\ d\leq
y\end{array}\right.$ (5.3)

and

$h(x,y)=\left\\{\begin{array}[]{ll}\dfrac{p+qd^{2}f_{XY}(x,y)}{p(d-x)(d-y)+qd^{2}R_{2}(x,y)},&\quad
0\leq x\leq d,\ 0\leq y\leq d\\\ \dfrac{f_{XY}(x,y)}{R_{2}(x,y)},&\quad d\leq
x\ \text{or}\ d\leq y\end{array}\right.$ (5.4)

Graphical plots are important to identify whether the models useful for
specific data sets for which empirical plots are available. For the
instantaneous failure case, the graph of the MBW density function, survival
function and hazard rate function are given for specific parameter values in
the following way.

  * 1.

Density function: The plot of the density functions for the parameter values
$\alpha_{1}=2$, $\beta_{1}=1$, $\alpha_{2}=2$, $\beta_{2}=1$, $a=1$, $b=1$,
$\rho=0.5$, $d=0.4$ and $p=0.2,\ 0.5,\ 0.8$ are given in Figure 1.

  * 2.

Survival function: Corresponding to the density functions given in Figure 1,
the plot of the survival functions are given in Figure 2.

  * 3.

Hazard function: The plot of the hazard function are given in the Figures 3,
4, 5 and 6 for different parameters. We can observe that the graph is similar
to the Weibull hazard function outside the rectangle $[0,d]\times[0,d]$.
Inside the rectangle, the behavior of the hazard function has more than one
direction where it takes increasing, decreasing, and bathtub shape.

Figure 1: Plot of MBW Density function with parameters: $\alpha_{1}=2$,
$\beta_{1}=1$, $\alpha_{2}=2$, $\beta_{2}=1$, $a=1$, $b=1$, $\rho=0.5$,
$d=0.4$, $p=0.2,\ 0.5,\ 0.8$. Figure 2: Plot of survival function with
parameters: $\alpha_{1}=2$, $\beta_{1}=1$, $\alpha_{2}=2$, $\beta_{2}=1$,
$a=1$, $b=1$, $\rho=0.5$, $d=0.4$, $p=0.2,\ 0.5,\ 0.8$. Figure 3: Plot of
hazard function with parameters: $\alpha_{1}=0.5$, $\beta_{1}=1$,
$\alpha_{2}=0.5$, $\beta_{2}=1$, $a=1$, $b=1$, $\rho=0.5$, $d=0.4$, $p=0.3$.
Figure 4: Plot of hazard function with parameters: $\alpha_{1}=0.5$,
$\beta_{1}=1$, $\alpha_{2}=0.5$, $\beta_{2}=1$, $a=1$, $b=1$, $\rho=0.5$,
$d=0.4$, $p=0.08$. Figure 5: Plot of hazard function with parameters:
$\alpha_{1}=3$, $\beta_{1}=1$, $\alpha_{2}=3$, $\beta_{2}=1$, $a=1$, $b=1$,
$\rho=0.5$, $d=0.4$, $p=0.08$. Figure 6: Plot of hazard function with
parameters: $\alpha_{1}=3$, $\beta_{1}=1$, $\alpha_{2}=3$, $\beta_{2}=1$,
$a=1$, $b=1$, $\rho=0.5$, $d=0.4$, $p=0.5$.

##  6 Estimation parameters

We estimate the unknown parameters of the proposed model using the MLE
combined with the DBSCAN algorithm as follows. The joint PDF of the modified
bivariate Weibull distribution allowing early and instantaneous failures is
given by

$f(x,y)=\left\\{\begin{array}[]{ll}\dfrac{p}{d^{2}}+qf_{XY}(x,y),&\quad 0\leq
x\leq d,\ 0\leq y\leq d\\\
qf_{XY}(x,y),&\quad\text{Otherwise,}\end{array}\right.$ (6.1)

where $f_{XY}$ is the joint PDF of bivariate Weibull distribution given in
(3.4). The unknown parameters of the model are given by
$\Theta=\\{\alpha_{1},\beta_{1},\alpha_{2},\beta_{2},\bm{\rho},d,p\\}$ with
$\bm{\rho}$ representing copula parameters. To estimate the parameters
$\Theta$, we combine the MLE technique and the DBSCAN clustering method in the
following way.
Let $(x_{i},y_{i}),\ i=1,2,\cdots,n$ be the $n$ sample points from the
population with PDF given by (6.1). Then, the likelihood function is

$L(\Theta|\textbf{x,y})=\prod_{0\leq x_{i}\leq d,\ 0\leq y_{i}\leq
d}\left[\dfrac{p}{d^{2}}+qf_{XY}(x_{i},y_{i})\right]\times\prod_{x_{i}>d\ or\
y_{i}>d}qf_{XY}(x_{i},y_{i}).$ (6.2)

and the corresponding log-likelihood function is

$\begin{split}l(\Theta|\textbf{x,y})&=log[L(\Theta|\textbf{x,y})]\\\
&=\sum_{0\leq x_{i}\leq d,\ 0\leq y_{i}\leq
d}log\left[\dfrac{p}{d^{2}}+qf_{XY}(x_{i},y_{i})\right]+\sum_{x_{i}>d\ or\
y_{i}>d}log\left[qf_{XY}(x_{i},y_{i})\right].\end{split}$ (6.3)

We usually maximize the log-likelihood function to find the estimates of
unknown parameters of the models by differentiating the likelihood or log-
likelihood function with respect to the parameters and solving the equations
simultaneously. However, we do not get a definite solution for d by the said
approach due to the particular nature of the likelihood function concerning
the corresponding variable. Now, we see that the log-likelihood function
decreases as the value of d increases. So, the maximum log-likelihood function
value occurred at the minimum possible value of d. Now, the contribution of d
in the likelihood function is non-zero when d satisfies the following
inequalities: $0<x_{i}<d$ and $0<y_{i}<d$ where $x_{i}$ and $y_{i}$ are
independently generated from uniform distribution in [0,d]. So, we classify
the data using the DBSCAN clustering method given by Ester et al., (1996). The
DBSCAN parameters are estimated using the methodologies provided in Schubert
et al., (2017) and Rahmah and Sitanggang, (2016). Then two clusters are
formed; one is C1, generated from the bivariate uniform distribution, and
another C2 from bivariate Weibull distribution. Now, we find the maximum
likelihood estimate of d given by,

$\hat{d}=\max\\{\max\\{x_{i},y_{i}\\}|(x_{i},y_{i})\in C1\\}$ (6.4)

Then, we use the usual approach of finding the MLEs of other parameters by
maximizing the log-likelihood function using the estimated value $\hat{d}$. In
the numerical examples given in Section 8, we use R-packages “optim” and
“dbscan” for estimating the model parameters using the proposed method.

##  7 Algorithm for Data Simulation

In this section, we discuss an algorithm for simulating data (x,y) from the
MBW distribution given in (4.1). The steps are as follows:

  1. 1.

Generates a value z from binomial distribution B(1,p).

  2. 2.

If $z=1$, then Generate x and y independently from the uniform distribution
$U(0,d)$.

  3. 3.

If $z=0$, then generates (x,y) from the bivariate Weibull distribution given
in (3.4), in the following way;

    1. (a)

Generate $u$ and $t$ independently from uniform $(0,1)$.

    2. (b)

Let $v=C^{-1}_{u}(t)$, where $C^{-1}_{u}$ is the quasi-inverse of $C_{u}$.
Here, $C_{u}$ is defined as, $C_{u}:=P(V\leq v|U=u)=\dfrac{\partial
C(u,v)}{\partial u}$, where C is a bivariate copula.
In the case of the GFGM copula, we can find $v$ by solving the following
equation using the Newton Raphson iteration method.

$v+\rho bu^{b-1}v^{b}(1-u)^{a}(1-v)^{a}-\rho
au^{b}v^{b}(1-u)^{a-1}(1-v)^{a}-t=0$

    3. (c)

Set $x=F_{1}^{-1}(u)$ and $y=F_{2}^{-1}(v)$, where $F_{1}$ and $F_{2}$ are the
CDF of univariate Weibull distribution with parameters
$(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ respectively.

To simulate n-data points using this algorithm we repeat these three steps n
times.

##  8 Numerical Example

Here, we provide numerical examples using simulated data to illustrate and
test the proposed methodologies and apply them to real data.

###  8.1 Simulation Studies

We simulate data $(x_{i},y_{i})$, for $i=1,2,3,\cdots n$, using the algorithm
proposed in Section 7 from modified bivariate Weibull distribution (4.1) with
parameter values assumed as $\Theta=(\alpha_{1}=4,\
\beta_{1}=1.5,\alpha_{2}=3.5,\ \beta_{2}=5,\ \rho=0.6,\ d=0.1,\ p=0.3)$. Three
different samples sizes $n=100,\ 200,$ and $300$ are considered for our
simulation studies. We fit the simulated data using the proposed modified
bivariate Weibull distribution (4.1) assuming the parameter vector $\Theta$ as
unknown.

The unknown parameter vector
$\Theta=\\{\alpha_{1},\beta_{1},\alpha_{2},\beta_{2},\rho,d,p\\}$ is estimated
using the proposed methodologies given in Section 6. We first find the
estimate $\hat{d}$ of $d$ using the DBSCAN clustering algorithm and the
statistic given in (6.4). For the sample sizes $n=100,200$, and $300$, the
estimated parameter values of DBSCAN $[Minpts=4,eps=0.45]$,
$[Minpts=4,eps=0.35]$, and $[Minpts=4,eps=0.25]$ respectively. Then, the log-
likelihood function is maximized to find the MLEs of the other model
parameters using the value of $d$ as $\hat{d}$. For each sample size
$n\in\\{100,200,300\\}$, we simulate 2000 data sets and estimate the model
parameters by the above methods. Then, we find the sample mean, Bias, Mean
Square Error (MSE), bootstrap standard error (BSE), and 95% bootstrap
confidence interval (BCI) using 2000 MLEs of each parameter for comparing the
estimated values with the true parameter values of the model. A 95% confidence
interval is also computed using the asymptotic results of MLE, and the
corresponding coverage probabilities are also found. The coverage probability
is considered as a percentage of times that the confidence interval includes
the respective true parameter values. To calculate the confidence interval of
$d$ for each sample data, we consider the univariate data
$\\{x_{i}|(x_{i},y_{i})\in C_{1}\\}\cup\\{y_{i}|(x_{i},y_{i})\in C_{1}\\}$ of
size $2n$ which follows uniform distribution on $[0,d]$, where $C_{1}$ is the
cluster near origin. The MSE and the Bias are calculated using the formula

$\begin{split}Bias(\hat{\theta})=\dfrac{1}{N}\sum_{i=1}^{N}\hat{\theta}^{(i)}-\theta\\\
MSE(\hat{\theta})=\dfrac{1}{N}\sum_{i=1}^{N}\left(\hat{\theta}^{(i)}-\theta\right)^{2}\end{split}$
(8.1)

where $\hat{\theta}$ denotes the estimator of the parameter $\theta$ and
$\hat{\theta}^{(i)}$ denotes the estimate of $\theta$ for $i$th data set,
where $N$ is the number of data generated with sample size $n$.

The sample mean, coverage probability, Bootstrap confidence interval (BCI),
MSE and Bias using 2000 estimated parameter values with sample sizes
$n=100,200,300$ are given in Table 1, Table 2 and Table 3 respectively. We see
that the Bias and MSE are small, i.e., the estimated parameter values by the
proposed methodologies are close to the true parameter values, and the
coverage probabilities are approximately the same as the confidence
coefficients for each scenario with varying sample sizes. The bootstrap
confidence intervals provide expected results for estimating model parameters
using the proposed methodologies.

Parameters | Sample Mean | CP | BSE | BCI | MSE | Bias
---|---|---|---|---|---|---
$\alpha_{1}=4$ | 4.0823 | 0.95 | 0.4360 | (3.4055 4.9091) | 0.1967 | 0.0823
$\beta_{1}=1.5$ | 1.4982 | 0.94 | 0.0581 | (1.4070 1.5939) | 0.0033 | -0.0017
$\alpha_{2}=3.5$ | 3.5628 | 0.94 | 0.3936 | (2.9375 4.3292) | 0.1588 | 0.0628
$\beta_{2}=5$ | 4.9818 | 0.94 | 0.2711 | (4.6446 5.3622) | 0.0738 | -0.0181
$\rho=0.6$ | 0.5970 | 0.93 | 0.0789 | (0.4342 0.7394) | 0.0062 | -0.0029
$d=0.1$ | 0.1005 | 0.94 | 0.0284 | (0.0934 0.0999) | 0.0008 | 0.0005
$p=0.3$ | 0.2969 | 0.95 | 0.0507 | (0.2099 0.3899) | 0.0025 | -0.0030
Table 1: Sample mean of estimated parameters, coverage probability(CP), 95% bootstrap confidence interval(BCI), MSE and Bias for 100 data points Parameters | Sample Mean | CP | BSE | BCI | MSE | Bias
---|---|---|---|---|---|---
$\alpha_{1}=4$ | 4.0068 | 0.94 | 0.4399 | (3.4956 4.6019) | 0.1935 | 0.0068
$\beta_{1}=1.5$ | 1.4935 | 0.94 | 0.0852 | (1.4277 1.5654) | 0.0072 | -0.0064
$\alpha_{2}=3.5$ | 3.5042 | 0.94 | 0.3918 | (3.0669 4.0266) | 0.1534 | 0.0042
$\beta_{2}=5$ | 4.9535 | 0.94 | 0.4224 | (4.7157 5.2445) | 0.1804 | -0.0464
$\rho=0.6$ | 0.5981 | 0.93 | 0.0679 | (0.4851 0.6995) | 0.0046 | -0.0018
$d=0.1$ | 0.1033 | 0.92 | 0.0383 | (0.0969 0.0999) | 0.0014 | 0.0033
$p=0.3$ | 0.2988 | 0.93 | 0.0516 | (0.2349 0.3700) | 0.0026 | -0.0011
Table 2: Sample mean of estimated parameters, coverage probability(CP), 95% bootstrap confidence interval(BCI), MSE and Bias for 200 data points Parameters | Sample Mean | CP | BSE | BCI | MSE | Bias
---|---|---|---|---|---|---
$\alpha_{1}=4$ | 4.0225 | 0.93 | 0.4264 | (3.3936 4.7940) | 0.1822 | 0.0225
$\beta_{1}=1.5$ | 1.4939 | 0.94 | 0.0660 | (1.4055 1.5787) | 0.0043 | -0.0060
$\alpha_{2}=3.5$ | 3.5349 | 0.94 | 0.3782 | (2.9488 4.2028) | 0.1442 | 0.0349
$\beta_{2}=5$ | 4.9768 | 0.94 | 0.2719 | (4.6402 5.3101) | 0.0744 | -0.0231
$\rho=0.6$ | 0.6010 | 0.94 | 0.0722 | (0.4573 0.7296) | 0.0052 | 0.0010
$d=0.1$ | 0.1025 | 0.95 | 0.0360 | (0.0946 0.0999) | 0.0013 | 0.0025
$p=0.3$ | 0.2984 | 0.94 | 0.0451 | (0.2199 0.3800) | 0.0020 | -0.0015
Table 3: Sample mean of estimated parameters, coverage probability(CP), 95%
bootstrap confidence interval(BCI), MSE and Bias for 300 data points

###  8.2 Application to Vannman wood data

An experiment was done by Vannman (1991) to compare two chemical processes for
drying wooden boards. Two different chemical processes are used to dry a batch
of wooden boards under the same climate conditions. In the experiment, two
chemical processes are compared, and damage to the board is measured as the
percentage of the checking area given in Vannman, (1991). As not all the
boards are checked, the sample contains some zero observation. This
observation corresponds to instantaneous failure or early failure. The
reproduced data set of the experiment on two batches of 36 boards by using two
different schedules are given in Table 4.

Sl. no. | schedule 1 $(y_{1})$ | schedule 2 $(y_{2})$ | Sl. no | schedule 1 $(y_{1})$ | schedule 2 $(y_{2})$
---|---|---|---|---|---
1 | 0 | 0 | 19 | 0.82 | 0,02
2 | 0 | 0 | 20 | 1.15 | 0.02
3 | 0 | 0 | 21 | 1.23 | 0.04
4 | 0 | 0 | 22 | 1.40 | 0.09
5 | 0 | 0 | 23 | 3.00 | 0.23
6 | 0 | 0 | 24 | 3.23 | 0.26
7 | 0 | 0 | 25 | 4.03 | 0.37
8 | 0 | 0 | 26 | 4.20 | 0.93
9 | 0 | 0 | 27 | 5.04 | 0.94
10 | 0 | 0 | 28 | 5.36 | 1.02
11 | 0 | 0 | 29 | 6.12 | 2.23
12 | 0 | 0 | 30 | 6.79 | 2.79
13 | 0 | 0 | 31 | 7.90 | 3.93
14 | 0.08 | 0 | 32 | 8.27 | 4.47
15 | 0.32 | 0 | 33 | 8.62 | 5.12
16 | 0.38 | 0 | 34 | 9.50 | 5.19
17 | 0.46 | 0 | 35 | 10.15 | 5.39
18 | 0.71 | 0.02 | 36 | 10.58 | 6.83
Table 4: Vannman experiment data(1991)

We consider three models, M1, M2, and M3, to fit the data for comparing the
proposed model M3 with existing other models M1 and M2, frequently used in the
literature. Model M1 assumes independent bivariate responses, and Model M2
considers dependent structure using a copula, ignoring the presence of
instantaneous and early failure in the data. As the data contains many $(0,0)$
observations, it will be realistic to fix the shape parameter value $\alpha=1$
in the two parameters Weibull distribution, i.e., exponential distribution,
otherwise for $\alpha\neq 1$, the log-likelihood value is not finite for model
M1 and M2. For the same reason, Gaussian copula, t-copula, and Archimedean
copulas are not applicable for the dependence structure in model M2, so we
used the FGM copula with the parameter values $a=1,\ b=1$. However, the
proposed model M3 can take two-parameter Weibull distributions as their
marginals and any copulas for the dependence structure. Here, we fit the data
using Gaussian, t-copula, Archimedean copula, and FGM copula, and the minimum
AIC value selects the FGM copula for Model M3.

The unknown parameters of Models M1 and M2 are estimated using the MLE, and
the proposed methodologies given in Section 6 are applied to estimate the
unknown parameters of Model M3. The parameter $d$ in Model M3 is found using
the estimated DBSCAN parameter values as $Minpts=4$, and $eps=1.6$. The
estimates, likelihood values, and Akaike Information Criterion (AIC) values
for Models M1, M2, and M3 are given in Table 5, Table 6, and Table 7. The AIC
values of Models M1, M2, and M3 are 228.4698, 195.8067, and 165.2824,
respectively, indicating that Model M3 is the best-fitted model to the data
with the minimum AIC values as a criterion. We also observed that Model M3
provides a significantly better fit than Models M1 and M2 by testing the
hypotheses using the deviance statistic, rejecting the Null hypothesis with
P-values less than 0.0001. Hence, researchers may find a better performance by
using the proposed model to fit bivariate continuous responses in the presence
of instantaneous and early failure in the data.

Parameters | Estimates | SE | p-value
---|---|---|---
$\beta_{1}$ | 2.759445 | 0.4599071 | $<0.0001$
$\beta_{2}$ | 1.108056 | 0.1846755 | $<0.0001$
Likelihood value | -112.2349 |  |
AIC | 228.4698 |  |
Table 5:  Parameter estimates, likelihood values, and AIC from model M1 for real data  Parameters | Estimates | SE | p-value
---|---|---|---
$\beta_{1}$ | 2.8934625 | 0.4586 | $<0.0001$
$\beta_{2}$ | 1.0333664 | 0.1642 | $<0.0001$
$\rho$ | 0.9994815 | 0.00000002 | $<0.0001$
Likelihood value | -94.90337 |  |
AIC | 195.8067 |  |
Table 6: Parameter estimates, likelihood values, and AIC from model M2 for real data  Parameters | Estimates | SE | p-value
---|---|---|---
$\alpha_{1}$ | 2.691784 | 0.8092 | $<0.0001$
$\beta_{1}$ | 7.739789 | 0.1405 | $<0.0001$
$\alpha_{2}$ | 1.000698 | 0.3492 | $<0.0001$
$\beta_{2}$ | 3.309285 | 0.00000002 | $<0.0001$
$\rho$ | 0.9849288 | 0.0164 | $<0.0001$
$d$ | 1.4 | 0.0017 | $<0.0001$
$p$ | 0.5784499 | 0.0847 | $<0.0001$
Likelihood value | -75.6412 |  |
AIC | 165.2824 |  |
Table 7: Parameter estimates, likelihood values, and AIC from model M3 for
real data

##  9 Conclusion

We proposed a modified bivariate Weibull distribution for modeling data
allowing early and instantaneous failure observations. The bivariate Weibull
distribution is obtained using a copula, assuming the marginals are
distributed as two parameters Weibull distribution. The survival and hazard
functions are provided, and the explicit forms of the functions are given in
the particular case of GFGM copula. We used MLE combined with DBSCAN
clustering algorithm to estimate the unknown model parameters. Numerical
examples are shown to illustrate and test the proposed method using simulated
data. The coverage probability, bootstrap standard error, bootstrap confidence
intervals, mean square error, and bias are calculated to test the performance
of the proposed methodologies. The proposed model is applied to real data and
compared with existing models in the literature. We observed that the modified
Weibull distribution outperforms with respect to the AIC as a criterion and
provides a significantly better fit by the chi-square test than other existing
models in the literature.

We have used the MLE combined with the DBSCAN clustering algorithm to estimate
the parameters. Researchers may use different estimation methods, such as the
Expectation-Maximization algorithm, Bayesian techniques, etc., to estimate the
model’s unknown parameters. One of the future directions may be to apply and
compare different estimation techniques to have an optimum performance of the
model.

## References

  * Almetwally et al., (2020) Almetwally, E. M., Muhammed, H. Z., and El-Sherpieny, E.-S. A. (2020).  Bivariate weibull distribution: properties and different methods of estimation.  Annals of Data Science, 7(1):163–193.
  * Arbenz, (2013) Arbenz, P. (2013).  Bayesian copulae distributions, with application to operational risk management—some comments.  Methodology and computing in applied probability, 15(1):105–108.
  * Cheng and Sheu, (2016) Cheng, Y.-F. and Sheu, S.-H. (2016).  Robust estimation for weibull distribution in partially accelerated life tests with early failures.  Quality and Reliability Engineering International, 32(7):2207–2216.
  * El-Gohary and El-Morshedy, (2015) El-Gohary, A. and El-Morshedy, M. (2015).  Bivariate exponentiated modified weibull extension.  arXiv preprint arXiv:1501.03528.
  * Ester et al., (1996) Ester, M., Kriegel, H.-P., Sander, J., Xu, X., et al. (1996).  A density-based algorithm for discovering clusters in large spatial databases with noise.  In kdd, volume 96, pages 226–231.
  * Joe, (2014) Joe, H. (2014).  Dependence modeling with copulas.  CRC press, Florida.
  * Lai et al., (2007) Lai, C., Khoo, M. B., Muralidharan, K., and Xie, M. (2007).  Weibull model allowing nearly instantaneous failures.  Journal of Applied Mathematics and Decision Sciences, 2007.
  * Lai and Xie, (2000) Lai, C. and Xie, M. (2000).  A new family of positive quadrant dependent bivariate distributions.  Statistics & probability letters, 46(4):359–364.
  * Muralidharan and Bavagosai, (2018) Muralidharan, K. and Bavagosai, P. (2018).  A new weibull model with inliers at zero and one based on type-ii censored samples.  Journal of the Indian Society for Probability and Statistics, 19(1):121–151.
  * Muralidharan and Bavagosai, (2021) Muralidharan, K. and Bavagosai, P. (2021).  Some inferences on tests for inliers.  International Journal of Statistics and Reliability Engineering, 8(2):243–263.
  * Muralidharan and Khabia, (2011) Muralidharan, K. and Khabia, A. (2011).  A modified pareto distribution.  J Indian Stat Assoc, 49:73–90.
  * Muralidharan and Lathika, (2006) Muralidharan, K. and Lathika, P. (2006).  Analysis of instantaneous and early failures in weibull distribution.  Metrika, 64(3):305–316.
  * Muralidharan et al., (2011) Muralidharan, K., Parikh, R., and Lai, C. (2011).  A bayesian analysis on weibull model allowing nearly instantaneous failures.  Reliability: Theory & Applications, 6(4 (23)):140–148.
  * Murthy et al., (2004) Murthy, D. P., Xie, M., and Jiang, R. (2004).  Weibull models, volume 505.  John Wiley & Sons.
  * Nelsen, (2007) Nelsen, R. B. (2007).  An introduction to copulas.  Springer Science & Business Media.
  * Peres et al., (2018) Peres, M. V. d. O., Achcar, J. A., and Martinez, E. Z. (2018).  Bivariate modified weibull distribution derived from farlie-gumbel-morgenstern copula: a simulation study.  Electronic Journal of Applied Statistical Analysis, 11(2):463–488.
  * Pham and Lai, (2007) Pham, H. and Lai, C.-D. (2007).  On recent generalizations of the weibull distribution.  IEEE transactions on reliability, 56(3):454–458.
  * Rahmah and Sitanggang, (2016) Rahmah, N. and Sitanggang, I. S. (2016).  Determination of optimal epsilon (eps) value on dbscan algorithm to clustering data on peatland hotspots in sumatra.  In IOP conference series: earth and environmental science, volume 31, page 012012. IOP Publishing.
  * Schubert et al., (2017) Schubert, E., Sander, J., Ester, M., Kriegel, H. P., and Xu, X. (2017).  Dbscan revisited, revisited: why and how you should (still) use dbscan.  ACM Transactions on Database Systems (TODS), 42(3):1–21.
  * Sklar, (1959) Sklar, M. (1959).  Fonctions de repartition an dimensions et leurs marges.  Publ. inst. statist. univ. Paris, 8:229–231.
  * Song and Song, (2007) Song, X.-K. and Song, P. X.-K. (2007).  Correlated data analysis: modeling, analytics, and applications.  Springer Science & Business Media, New York.
  * Vannman, (1991) Vannman, K. (1991).  Comparing samples from nonstandard mixtures of distributions with applications to quality comparison of wood. research report 1991:2.  submitted to division of quality technology, Lulea University, Lulea, Sweden.
  * Withana Gamage et al., (2020) Withana Gamage, P. W., Chaudari, M., McMahan, C. S., Kim, E. H., and Kosorok, M. R. (2020).  An extended proportional hazards model for interval-censored data subject to instantaneous failures.  Lifetime data analysis, 26(1):158–182.


# QQQ: Quality Quattuor-Bit Quantization for Large Language Models

Ying Zhang1, Peng Zhang1, Mincong Huang1, Jingyang Xiang1,2, Yujie Wang1
Chao Wang1, Yineng Zhang1, Lei Yu1, Chuan Liu, Wei Lin
1Meituan 2Zhejiang University
<EMAIL_ADDRESS>
###### Abstract

Quantization is a proven effective method for compressing large language
models. Although popular techniques like W8A8 and W4A16 effectively maintain
model performance, they often fail to concurrently speed up the prefill and
decoding stages of inference. W4A8 is a promising strategy to accelerate both
of them while usually leads to a significant performance degradation. To
address these issues, we present QQQ, a Quality Quattuor-bit Quantization
method with 4-bit weights and 8-bit activations. QQQ employs adaptive
smoothing and Hessian-based compensation, significantly enhancing the
performance of quantized models without extensive training. Furthermore, we
meticulously engineer W4A8 GEMM kernels to increase inference speed. Our
specialized per-channel W4A8 GEMM and per-group W4A8 GEMM achieve impressive
speed increases of 3.67$\times$ and 3.29 $\times$ over FP16 GEMM. Our
extensive experiments show that QQQ achieves performance on par with existing
state-of-the-art LLM quantization methods while significantly accelerating
inference, achieving speed boosts up to 2.24 $\times$, 2.10$\times$, and
1.25$\times$ compared to FP16, W8A8, and W4A16, respectively. Code is
available at https://github.com/HandH1998/QQQ.

QQQ: Quality Quattuor-Bit Quantization for Large Language Models

Ying Zhang1, Peng Zhang1, Mincong Huang1, Jingyang Xiang1,2, Yujie Wang1 Chao
Wang1, Yineng Zhang1, Lei Yu1, Chuan Liu, Wei Lin 1Meituan 2Zhejiang
University<EMAIL_ADDRESS>

##  1 Introduction

In recent years, large language models (LLMs) have been widely applied in
various tasks and shown excellent performance, such as OPT (Zhang et al.,
2022) and LLaMA series Touvron et al. (2023a, b). The LLMs usually have large
number of parameters and long inference time, making it inapplicable to
resource-limited devices and real-time scenarios. Therefore, it is crucial to
reduce LLMs’ storage and computation overhead while retaining their
performance.

Quantization is a crucial technique for reducing the memory and computational
requirements of LLMs. It generally involves two main strategies: weight-
activation and weight-only quantization. The W8A8 weight-activation
quantization method, as demonstrated by SmoothQuant (Xiao et al., 2023),
compresses both weights and activations to 8-bit, leveraging specialized
hardware like INT8 Tensor Cores to increase computational throughput. This
method excels in enhancing the compute-bound prefill phase of LLM inference.
Weight-only quantization, such as the W4A16 method seen in GPTQ (Frantar et
al., 2022) and AWQ (Lin et al., 2023), reduces the weight precision to 4-bit
while maintaining 16-bit activations. This approach effectively cuts memory
bandwidth requirements and is particularly beneficial during the memory-bound
decoding phase of LLM inference. Although these methods successfully retain
model performance, they struggle to simultaneously accelerate both the prefill
and decoding stages of LLM inference. The W4A8 weight-activation quantization
approach, which can expedite both inference phases by utilizing efficient
hardware and reducing memory bandwidth, often leads to a significant decrease
in model performance. Therefore, the pursuit of a viable W4A8 quantization
method that balances performance with speed remains an active area of
research.

In response to these challenges, we introduce QQQ, a quality W4A8 quantization
approach that excels in both model performance and inference speed. QQQ
employs an adaptive smoothing mechanism for activation channels with
significant outliers, while leaving other channels intact. This targeted
smoothing allows for more effective activation quantization. Following the
smoothing process, QQQ utilizes a Hessian-based compensation technique to
offset the loss from weight quantization. Through these innovative strategies,
QQQ significantly improves the performance of quantized models. To achieve
fast inference speed, we meticulously design specialized W4A8 GEMM kernels.
These kernels are specifically tailored for both per-channel and per-group
weight quantization, drastically enhancing throughput across various batch
sizes. Our contributions compared to previous works are threefold:

  * •

We present QQQ, a quality W4A8 quantization method, that delivers performance
on par with existing state-of-the-art LLM quantization methods.

  * •

We develop innovative W4A8 GEMM kernels tailored for both per-channel weight
and per-group weight quantization. These specialized kernels significantly
enhance speed, with per-channel W4A8 GEMM and per-group W4A8 GEMM achieving
speedups of 3.67$\times$ and 3.29$\times$ over FP16 GEMM, respectively.

  * •

We integrate QQQ into vLLM (Kwon et al., 2023), achieving impressive speed
boosts of 2.24$\times$, 2.10$\times$, and 1.25$\times$ over FP16, W8A8, and
W4A16 Marlin (Frantar and Alistarh, 2024) implementation, respectively.

##  2 Related Works

###  2.1 Large Language Models

In recent times, large language models (LLMs) such as GPT-3 (Brown et al.,
2020), GLM (Du et al., 2021), OPT (Zhang et al., 2022), and the LLaMA series
Touvron et al. (2023a, b) have revolutionized the field of natural language
processing by delivering unprecedented performance enhancements across various
tasks. Despite their impressive capabilities, LLMs come with a significant
computational footprint due to their billions of parameters, which takes long
inference time. This makes them less suitable for deployment on devices with
limited resources or in applications requiring real-time responses. For
example, the GPT-3 (Brown et al., 2020) model, with its colossal count of 175
billion parameters, demands a staggering 350 GB of memory for loading its
parameters in FP16 format. This necessitates the use of at least five A100-80G
GPUs solely for inference purposes.

###  2.2 Quantization

Current strategies for quantizing LLMs fall into two main categories: weight-
activation and weight-only quantization.

Weight-activation quantization. Weight-activation quantization compresses both
weights and activations to lower bit representations such as W8A8, W4A8, W6A6,
and W4A4. SmoothQuant (Xiao et al., 2023), LLM.int8() (Dettmers et al., 2022),
QUIK (Ashkboos et al., 2023) and Outlier Suppression+ (Wei et al., 2023) all
enhance quantized model performance by managing activation outliers.
LLM.int8() (Dettmers et al., 2022) and QUIK (Ashkboos et al., 2023) employ
mixed-precision decomposition, while SmoothQuant (Xiao et al., 2023) and
Outlier Suppression+ (Wei et al., 2023) use channel-wise scaling. RPTQ (Yuan
et al., 2023) mitigates the impact of range differences between channels by
rearranging the channels and quantizing them in clusters. OmniQuant (Shao et
al., 2023) proposes learnable weight clipping and learnable equivalent
transformation. Atom (Zhao et al., 2023) attains high accuracy by applying a
novel mixed-precision and fine-grained quantization process. QuaRot (Ashkboos
et al., 2024) removes outliers by rotating the inputs of the model using
randomized Hadamard transformations.

Weight-only quantization. Weight-only quantization compresses model weights to
lower bit precision (e.g., 4-bit or 2-bit) while keeping activations at higher
precision (16-bit). GPTQ (Frantar et al., 2022) adopts second-order
information to minimize precision loss. SpQR (Dettmers et al., 2023) and AWQ
(Lin et al., 2023) prioritize the preservation of important weights to reduce
quantization errors, with SpQR applying mixed-precision and AWQ optimizing
per-channel scaling. QuIP (Chee et al., 2024) further quantizes weights to
2-bit based on the insight that model parameters should ideally be incoherent.

##  3 Methodology

In this section, we detail our quantization method, referred to as QQQ. As
depicted in Figure 1, QQQ is a two-stage weight-activation quantization
method. It integrates adaptive smoothing and Hessian-based quantization
compensation to enhance the performance of quantized models. QQQ adaptively
smooths activation channels with significant outliers while leaving other
channels intact, enabling better activation quantization. After smoothing the
activations, QQQ tackles the challenges of weight quantization by employing
Hessian-based compensation, which effectively minimizes the loss incurred
during the quantization process. To expedite practical inference, QQQ
incorporates highly efficient W4A8 GEMMs tailored for both per-channel and
per-group quantization. The quantization process will be further elaborated in
the following subsections.

Figure 1: Framework of QQQ.

###  3.1 Adaptive Smoothing

We employ channel-wise smoothing to refine problematic activation channels
with outliers, enhancing their suitability for quantization. As indicated by
research (Wei et al., 2023), smoothing only the activation channels with
significant outliers while leaving other channels intact can enhance the
performance of quantized models. Therefore, we further implement an adaptive
smoothing strategy that identifies and selectively smooths only those channels
that display significant outliers.

Specifically, we introduce an outlier threshold $\sigma$, which is determined
through a grid search within the range
$[0,\mbox{max}(\mbox{abs}(\mathbf{X}))]$. Activation channels that exhibit
outliers exceeding this threshold $\sigma$ are then selectively smoothed:

$\mathcal{T}=\\{\,t\,|\,\mbox{max}(\mbox{abs}(\mathbf{X}_{:,t}))<\sigma\\},$
(1)

where $\mathcal{T}$ denotes the selected channel index for $\mathbf{X}$.

To optimize the smoothing parameters, we aim to reduce the squared error
between the product of activations and weights both before and after
quantization:

$\arg\min_{\mathbf{s}}\|Q(\mathbf{X}_{:,\mathcal{T}}\oslash\mathbf{s})Q(\mathbf{W}_{\mathcal{T},:}\odot\mathbf{s})-\mathbf{X}_{:,\mathcal{T}}\mathbf{W}_{\mathcal{T},:}\|_{2}^{2},$
(2)

where $Q$ denotes the quantization function; $\mathbf{W}$ denotes the
unquantized weight matrix; $\mathbf{X}$ denotes the unquantized activation
matrix; $\mathbf{s}$ denotes the smoothing parameter; $\oslash$ and $\odot$
are element-wise division and multiplication, respectively.

###  3.2 Hessian-based Quantization Compensation

After transferring the quantization challenge from activations to weights
through adaptive smoothing, it becomes crucial to enhance the quantization
technique for weights. We adopt a layer-wise quantization framework GPTQ
(Frantar et al., 2022), which advocates for a precise and efficient Hessian-
based quantization compensation technique. This method involves an application
of two key formulas, Eq.(3) and Eq.(4), which are instrumental in the
iterative quantization of weights. Specifically, this algorithm adjusts the
subset $F$ of full-precision weights by an increment $\delta_{F}$,
compensating for the quantization error induced by the quantized weights
$Q(\textbf{W}_{i})$.

$\mathbf{W}_{i}=\mathop{\mbox{argmin}}\limits_{\mathbf{W}_{i}}\frac{(Q(\textbf{W}_{i})-\textbf{W}_{i})^{2}}{[\mathbf{H}_{F}^{-1}]_{ii}}$
(3)
$\bm{\delta}_{F}=-\frac{\textbf{W}_{i}-Q(\textbf{W}_{i})}{[\mathbf{H}_{F}^{-1}]_{ii}}\cdot(\mathbf{H}_{F}^{-1})_{:,i}$
(4)

where $\textbf{W}_{i}$ denotes the $i$-row of weight matrix $\mathbf{W}$,
$\mathbf{H}_{F}=2\textbf{X}_{F}^{\mathsf{T}}\textbf{X}_{F}$ denotes the
Hessian matrix.

###  3.3 Novel W4A8 GEMM

We carefully develop innovative W4A8 GEMMs tailored to optimize the inference
speed. It’s worth noting that our design exclusively supports symmetric
quantization due to its superior hardware efficiency compared to asymmetric
quantization.

The W4A8 GEMM design introduces two critical considerations beyond
conventional GEMM implementations. Firstly, as GPUs support GEMM operations
only for identical data types, it’s necessary to convert INT4 weights into
INT8 to utilize the capabilities of the INT8 Tensor Cores. Secondly, to avoid
potential overflow during INT8 GEMM, the final output is typically stored in
INT32 format, which subsequently requires a dequantization setp from INT32 to
FP16.

To further reduce memory access, we integrate both the type conversion and
dequantization process within the W4A8 GEMM routine. We have different W4A8
GEMM configurations for per-channel weight quantization and per-group weight
quantization, as they have different quantization schemes, which will be
detailed in the following subsections.

Figure 2: The W4A8 GEMM for per-channel weight quantization. $\mathbf{A}$
denotes the activation matrix. $\mathbf{W}$ denotes the weight matrix.
$\mathbf{s_{A}}$ denotes the per-token quantization scale of activations.
$\mathbf{s_{W}}$ denotes the per-channel quantization scale of weights. The
superscript represents the data type. The green sections indicate data, and
the orange sections indicate operations.

####  3.3.1 Per-channel Weight Quantization

The specialized W4A8 GEMM process for per-channel weight quantization is
depicted in Figure 2. Our method involves an initial conversion of INT4
weights to INT8 before executing the INT8 GEMM, followed by a dequantization
step that transforms the INT32 output of the INT8 GEMM into a FP16 result. To
ensure high efficiency in the type conversion, we utilize the FastINT4toINT8
(Li et al., 2023) technique. This process involves positioning the INT4 weight
into the upper 4 bits of an INT8 by multiplying by 16, essentially performing
a left shift by 4 bits. Once the GEMM computation is complete, we scale down
the result by a factor of 16 to obtain the accurate output. This scaling
operation is seamlessly integrated into the dequantization step, Dequant, by
adjusting the weight quantization scale offline by the same factor.

Figure 3: The W4A8 GEMM for per-group weight quantization.
$\mathbf{s_{W_{c}}}$ represents the per-channel quantization scale converting
from dequantized FP16 weights back to INT8 weights. $\mathbf{s_{W_{g}}^{*}}$
represents the fused per-group scale converting from INT4 weights to INT8
weights. Figure 4: Design of FusedDequantQuant. $\mathbf{s_{W_{g}}}$
represents the per-group quantization scale converting from FP16 weights to
INT4 weights. The orange sections indicate online operations, and the blue
sections indicate offline operations.

####  3.3.2 Per-group Weight Quantization.

The W4A8 GEMM tailored for per-group weight quantization is detailed in Figure
3. This approach involves a more complex set of operations prior to the INT8
GEMM compared to the per-channel weight quantization. Although the conversion
from INT4 to INT8 weights is still necessary, additional steps are required
due to the unique per-group quantization scales that are not directly
compatible with standard GEMM procedures.

Initially, INT4 weights are converted to FP16 using the FastINT4toFP16 (Kim et
al., 2022) method. Subsequently, each group’s weight is dequantized by
multiplying it by its respective scale. To re-quantize into INT8, the result
is divided by the per-channel quantization scale. These two steps are
elegantly merged into a single operation, FusedDequantQuant, by dividing the
per-group scale by the per-channel scale, as illustrated in Figure 4. The
final stage of the process involves converting the FP16 results of
FusedDequantQuant back into INT8 format to proceed with the INT8 GEMM. Similar
to the per-channel method, the INT32 output from the INT8 GEMM undergoes
dequantization to FP16.

Figure 5: Design of FastFP16toINT8.

The conventional instruction for converting FP16 to INT8 is significantly
slow, which hinders the overall GEMM performance. To tackle this issue, we
introduce a fast conversion technique, FastFP16toINT8, as depicted in Figure
5. The conversion begins by adjusting the FP16 values, adding $128$ to shift
them into the range of $[0.0,255.0]$, which corresponds to the
representational span of UINT8. This is followed by an addition of $1024$,
effectively aligning the 8 bits of UINT8 within the lower segment of FP16’s
mantissa. The final step involves extracting these lower 8 bits from FP16 and
applying an XOR operation with $0\textrm{x}80$ to achieve the desired INT8
format. In practice, the fusion of the per-group quantization scale
multiplication with the addition of the magic number $128+1024$ is executed
through a single FMA (Fused Multiply-Add) instruction. Furthermore, the
process leverages two expedient bit manipulation instructions, PRMT and XOR,
to finalize the conversion.

#Bits | Method | LLaMA-1 | LLaMA-2 | LLaMA-3
---|---|---|---|---
7B | 13B | 30B | 7B | 13B | 70B | 8B
FP16 | - | 5.68 | 5.09 | 4.10 | 5.47 | 4.88 | 3.32 | 6.14
W8A8 | SmoothQuant | 5.78 | 5.17 | 4.28 | 5.58 | 4.93 | 3.39 | 6.25
W4A16 | GPTQ-g128 | 5.83 | 5.20 | 4.22 | 5.63 | 4.99 | 3.43 | 6.56
AWQ-g128 | 5.78 | 5.19 | 4.21 | 5.60 | 4.97 | 3.41 | 6.54
W4A4 | Atom-g128 | 6.25 | 5.52 | 4.61 | 6.12 | 5.31 | 3.73 | 7.76
W4A8 | QoQ | 5.93 | 5.28 | 4.34 | 5.75 | 5.12 | 3.52 | 6.89
QoQ-g128 | 5.89 | 5.25 | 4.28 | 5.70 | 5.08 | 3.47 | 6.76
QQQ | 6.19 | 5.43 | 4.61 | 5.95 | 5.21 | 3.68 | 7.41
QQQ-g128 | 5.87 | 5.24 | 4.30 | 5.71 | 5.01 | 3.50 | 6.64
Table 1: WikiText2 perplexity with 2048 sequence length. The lower is the better. ’-g128’ denotes applying per-group quantization on the weights and the group size is 128.  LLaMA-2 | #Bits | Method | PIQA | ARC-e | ARC-c | HellaSwag | WinoGrande | Avg.
---|---|---|---|---|---|---|---|---
7B | FP16 | - | 79.05 | 74.58 | 46.25 | 76.05 | 68.98 | 68.98
W4A4 | Atom-g128 | 75.14 | 52.99 | 38.40 | 69.37 | 62.75 | 59.73
W4A8 | QoQ | 77.64 | 72.81 | 43.60 | 74.00 | 68.03 | 67.22
QoQ-g128 | 78.07 | 73.32 | 44.80 | 74.98 | 68.59 | 67.95
| QQQ | 77.42 | 69.15 | 42.15 | 73.54 | 65.98 | 65.65
| QQQ-g128 | 78.51 | 72.94 | 44.37 | 74.53 | 67.01 | 67.47
13B | FP16 | - | 80.52 | 77.44 | 49.06 | 79.38 | 72.22 | 71.72
W4A4 | Atom-g128 | 76.50 | 57.49 | 42.32 | 73.84 | 67.40 | 63.51
W4A8 | QoQ | 79.71 | 75.97 | 48.38 | 77.80 | 70.96 | 70.56
QoQ-g128 | 79.43 | 77.06 | 48.81 | 78.35 | 70.48 | 70.83
| QQQ | 79.43 | 74.75 | 48.12 | 77.27 | 70.32 | 69.98
| QQQ-g128 | 79.98 | 76.64 | 48.55 | 78.63 | 71.82 | 71.13
Table 2: Zero-shot accuracy on five common sense tasks with 2048 sequence
length. The higher is the better.

##  4 Experiments

###  4.1 Setups

In this section, we introduce the experimental configurations in models,
datasets, algorithm and inference system.

Models and Datasets. We conduct the experiments using the LLaMA series Touvron
et al. (2023a, b). We evaluate the model performance on language modeling and
zero-shot tasks. For the language modeling aspect, we utilize the WikiText2
(Merity et al., 2016) dataset to measure perplexity. For the zero-shot tasks,
we extend the evaluation to a variety of datasets, including PIQA (Bisk et
al., 2020), HellaSwag (Zellers et al., 2019), WinoGrande (Sakaguchi et al.,
2021), and ARC (Clark et al., 2018), utilizing the lm_eval (Gao et al., 2023)
toolkit for a comprehensive analysis. Additionally, we select a subset of 128
sequences from the PILE (Gao et al., 2020) dataset to facilitate calibration.

Algorithm. QQQ is specifically designed to focus on per-token quantization for
activations and supports both per-channel and per-group quantization for
weights. To optimize inference speed, QQQ exclusively uses symmetric
quantization for both activations and weights, thus avoiding the additional
computational overhead associated with asymmetric quantization.

Inference System. The inference system that we have develop to support QQQ is
based on vLLM111https://github.com/vllm-project/vllm (Kwon et al., 2023)
v0.4.1. The system’s capabilities are further enhanced by the implementation
of W4A8 GEMM kernels, which are based on the advanced
Marlin222https://github.com/IST-DASLab/marlin (Frantar and Alistarh, 2024)
kernels. All experimental procedures are conducted on NVIDIA A100 80G GPUs
under PyTorch 2.1.2 with CUDA 11.8.

###  4.2 Model Performance

We perform a comprehensive comparison of QQQ method against leading LLM
quantization methods, including SmoothQuant (Xiao et al., 2023), GPTQ (Frantar
et al., 2022), AWQ (Lin et al., 2023), Atom (Zhao et al., 2023), and the
recently introduced W4A8 method QoQ (Lin et al., 2024). By default, weights
employ per-channel quantization. GPTQ, AWQ, and QoQ utilize asymmetric
quantization for weights, whereas SmoothQuant, Atom, and QQQ opt for symmetric
quantization. In terms of activation quantization, Atom stands out by
implementing per-group quantization, while the remaining methods, including
QQQ, employ per-token quantization. Across all methods under consideration,
activations apply symmetric quantization.

Table 1 presents a detailed perplexity comparison on the WikiText2 dataset.
The per-group QQQ demonstrates competitive performance with with W8A8
SmoothQuant, W4A16 GPTQ, and W4A16 AWQ across various models. For example,
per-group QQQ only increases perplexity by up to 0.13 on LLaMA-2-13B compared
with them. When benchmarked against W4A4 Atom, both the per-channel and per-
group QQQ consistently outperform Atom on all tested models. For instance,
per-group QQQ achieves a significant reduction in perplexity by 0.41 on
LLaMA-2-13B when compared to Atom. In comparison with W4A8 QoQ, per-group QQQ
outperforms per-group QoQ on the majority of models, exemplified by a 0.12
decrease in perplexity on LLaMA-3-8B. However, it is observed that per-channel
QQQ falls slightly behind per-channel QoQ. This divergence is hypothesized to
be due to QoQ’s application of asymmetric quantization on weights, as opposed
to the symmetric quantization employed by QQQ.

The zero-shot accuracy evaluation across five common sense tasks is showed in
Table 2. The QQQ method exhibits a significant performance edge over Atom
across all tasks and model sizes. For instance, the per-group QQQ achieves a
4.79% higher accuracy than Atom on the HellaSwag task when utilizing the
LLaMA-2-13B model. In comparison with QoQ, QQQ demonstrates comparable
results, especially when per-group quantization is applied. For instance, QQQ
shows a marginal improvement of 0.3% over QoQ with the LLaMA-2-13B model on
average.

Figure 6: Same-batch throughput comparison of quantized LLaMA-2 models under
various batch size.

###  4.3 Speedup

In the evaluation of inference efficiency, QQQ is benchmarked against a suite
of quantization methods including SmoothQuant333https://github.com/vllm-
project/vllm/pull/1508 (Xiao et al., 2023), AWQ (Lin et al., 2023), and GPTQ
(Frantar et al., 2022), all of which are integrated within the vLLM (Kwon et
al., 2023). For optimal inference speed, we leverage the Marlin (Frantar and
Alistarh, 2024) implementation of GPTQ. All methods are tested under
conditions that enable continuous batching and paged attention. The benchmarks
are conducted with an input sequence length of 1024 and an output sequence
length of 128.

The throughput comparison outlined in Figure 6 indicates that QQQ outperforms
the alternative methods across a range of batch sizes in most cases.
Specifically, QQQ attains a substantial speed increase, reaching up to
2.24$\times$, 2.10$\times$, 1.59$\times$, and 1.25$\times$ faster than FP16,
SmoothQuant, AWQ, and Marlin on the LLaMA-2-13B model, respectively. When
contrasted with W4A16 AWQ and W4A16 Marlin, QQQ displays a comparable
inference speed at smaller batch sizes but gains a significant advantage as
the batch size increases. Against W8A8 SmoothQuant, QQQ stands out with an
approximate 2.0$\times$ speedup for smaller batches and maintains a similar
level of acceleration even at a batch size of 64. Moreover, QQQ’s performance
enhancement is more pronounced on larger models, as evidenced by the greater
speedup achieved on LLaMA-2-13B compared to LLaMA-2-7B.

Figure 7: Speedup over PyTorch FP16 GEMM (Calling CUTLASS) of all GEMMs under
different numbers of input tokens. The weight matrix size is
($N=8192,K=21760$).

###  4.4 W4A8 GEMM Studies

We conduct a series of experiments to demonstrate the capabilities of our W4A8
GEMMs. By keeping the weight matrix size constant and varying the number of
input tokens, we evaluate the speed improvements over the PyTorch FP16 GEMM
benchmark. The results are depicted in Figure 7.

Our W4A8 GEMMs demonstrate impressive speedups, achieving up to 3.67$\times$,
2.51$\times$, and 1.71$\times$ faster performance compared to PyTorch FP16
GEMM, W8A8 GEMM, and W4A16 Marlin GEMM, respectively. Specifically, W4A8 GEMMs
show a significant advantage when the token count is below 64 compared to W8A8
GEMM. With an increasing number of tokens, the per-group W4A8 GEMM maintains a
similar level of efficiency, whereas the per-channel W4A8 GEMM tends to be
slower. This discrepancy is attributed to the higher type conversion overhead
in per-group W4A8 GEMM. Compared to Marlin GEMM, our GEMMs demonstrate
superior performance under the majority of token counts, especially per-group
GEMM. As the token count rises, Marlin GEMM’s performance gradually declines,
falling behind even the FP16 GEMM, while our GEMMs consistently outperform the
FP16 benchmark. This clearly indicates the superior efficiency of INT8 Tensor
Cores over FP16 Tensor Cores in compute-bound scenarios.

###  4.5 Ablation Studies

Granularity | Model | B | B+AS | B+AS+GPTQ
---|---|---|---|---
Per-channel | LLaMA-1-7B | 6.93 | 6.77 | 6.19
LLaMA-1-13B | 5.84 | 5.84 | 5.43
LLaMA-1-30B | 5.04 | 4.89 | 4.61
LLaMA-2-7B | 7.36 | 7.29 | 5.95
LLaMA-2-13B | 5.47 | 5.43 | 5.21
LLaMA-2-70B | 4.07 | 3.88 | 3.68
LLaMA-3-8B | 10.88 | 9.46 | 7.41
Per-group | LLaMA-1-7B | 6.16 | 5.96 | 5.87
LLaMA-1-13B | 5.43 | 5.32 | 5.24
LLaMA-1-30B | 4.51 | 4.34 | 4.30
LLaMA-2-7B | 5.86 | 5.77 | 5.71
LLaMA-2-13B | 5.10 | 5.06 | 5.01
LLaMA-2-70B | 3.67 | 3.52 | 3.50
LLaMA-3-8B | 7.11 | 6.96 | 6.64
Table 3: WikiText2 perplexity of LLaMA series when applying various
quantization techniques. ’B’ represents the Baseline, which refers to the
vanilla W4A8 method without any additional quantization techniques. ’AS’
denotes Adaptive Smoothing.

####  4.5.1 Impact of Quantization Techniques

Here we examine the impact of different quantization techniques on model
performance. We use the standard W4A8 method, which lacks any specialized
quantization techniques, as the baseline for comparison. The perplexity
results for the LLaMA series are summarized in Table 3. Overall, the
combination of adaptive smoothing with GPTQ proves to be the most effective
method, yielding the lowest perplexity scores. For instance, this approach
reduces perplexity by 0.47 in per-group quantization compared with the
baseline.

Figure 8: Speedup over PyTorch FP16 GEMM (Calling CUTLASS) of two per-group
W4A8 GEMMs under different numbers of input tokens. The weight matrix size is
($N=8192,K=21760$).

####  4.5.2 Impact of FastFP16toINT8

Here we explore the performance enhancements achieved in per-group W4A8 GEMM
through the implementation of the FastFP16toINT8 conversion technique. We
benchmark the optimized W4A8 GEMM against the conventional W4A8 GEMM, which
utilizes the standard FP16 to INT8 conversion instruction. As illustrated in
Figure 8, the GEMM utilizing FastFP16toINT8 demonstrates a significant
increase in speedup compared to its counterpart without this optimization,
across a range of input tokens. The FastFP16toINT8 technique can accelerate
per-group W4A8 GEMM by up to 2.02$\times$. This finding confirms that the FP16
to INT8 conversion process is a critical bottleneck, and our FastFP16toINT8
effectively enhances the efficiency of per-group W4A8 GEMM by mitigating this
constraint.

##  5 Conclusion

In this paper, we propose a quality W4A8 quantization method optimized for
hardware efficiency. QQQ leverages the power of adaptive smoothing combined
with Hessian-based quantization compensation to significantly improve the
performance of quantized models. In addition, we meticulously develop W4A8
GEMM kernels to maximize the inference speed. Through extensive experiments,
QQQ demonstrates superior performance over existing state-of-the-art LLM
quantization methods in terms of both performance and speed. We believe QQQ
can serve as a powerful tool that opens up access to the deployment and
application of LLMs.

## Limitations

QQQ method has two limitations: 1) The two-stage quantization method requires
more time for quantizing models. 2) The mixed-precision GEMMs only support
4-bit weights for now.

## References

  * Ashkboos et al. (2023) Saleh Ashkboos, Ilia Markov, Elias Frantar, Tingxuan Zhong, Xincheng Wang, Jie Ren, Torsten Hoefler, and Dan Alistarh. 2023.  Towards end-to-end 4-bit inference on generative large language models.  _arXiv preprint arXiv:2310.09259_.
  * Ashkboos et al. (2024) Saleh Ashkboos, Amirkeivan Mohtashami, Maximilian L Croci, Bo Li, Martin Jaggi, Dan Alistarh, Torsten Hoefler, and James Hensman. 2024.  Quarot: Outlier-free 4-bit inference in rotated llms.  _arXiv preprint arXiv:2404.00456_.
  * Bisk et al. (2020) Yonatan Bisk, Rowan Zellers, Jianfeng Gao, Yejin Choi, et al. 2020.  Piqa: Reasoning about physical commonsense in natural language.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 34, pages 7432–7439.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. 2020.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901.
  * Chee et al. (2024) Jerry Chee, Yaohui Cai, Volodymyr Kuleshov, and Christopher M De Sa. 2024.  Quip: 2-bit quantization of large language models with guarantees.  _Advances in Neural Information Processing Systems_ , 36.
  * Clark et al. (2018) Peter Clark, Isaac Cowhey, Oren Etzioni, Tushar Khot, Ashish Sabharwal, Carissa Schoenick, and Oyvind Tafjord. 2018.  Think you have solved question answering? try arc, the ai2 reasoning challenge.  _arXiv preprint arXiv:1803.05457_.
  * Dettmers et al. (2022) Tim Dettmers, Mike Lewis, Younes Belkada, and Luke Zettlemoyer. 2022.  Llm. int8 (): 8-bit matrix multiplication for transformers at scale, 2022.  _CoRR abs/2208.07339_.
  * Dettmers et al. (2023) Tim Dettmers, Ruslan Svirschevski, Vage Egiazarian, Denis Kuznedelev, Elias Frantar, Saleh Ashkboos, Alexander Borzunov, Torsten Hoefler, and Dan Alistarh. 2023.  Spqr: A sparse-quantized representation for near-lossless llm weight compression.  _arXiv preprint arXiv:2306.03078_.
  * Du et al. (2021) Zhengxiao Du, Yujie Qian, Xiao Liu, Ming Ding, Jiezhong Qiu, Zhilin Yang, and Jie Tang. 2021.  Glm: General language model pretraining with autoregressive blank infilling.  _arXiv preprint arXiv:2103.10360_.
  * Frantar and Alistarh (2024) Elias Frantar and Dan Alistarh. 2024.  Marlin: a fast 4-bit inference kernel for medium batchsizes.  https://github.com/IST-DASLab/marlin.
  * Frantar et al. (2022) Elias Frantar, Saleh Ashkboos, Torsten Hoefler, and Dan Alistarh. 2022.  Gptq: Accurate post-training quantization for generative pre-trained transformers.  _arXiv preprint arXiv:2210.17323_.
  * Gao et al. (2020) Leo Gao, Stella Biderman, Sid Black, Laurence Golding, Travis Hoppe, Charles Foster, Jason Phang, Horace He, Anish Thite, Noa Nabeshima, et al. 2020.  The pile: An 800gb dataset of diverse text for language modeling.  _arXiv preprint arXiv:2101.00027_.
  * Gao et al. (2023) Leo Gao, Jonathan Tow, Baber Abbasi, Stella Biderman, Sid Black, Anthony DiPofi, Charles Foster, Laurence Golding, Jeffrey Hsu, Alain Le Noac’h, Haonan Li, Kyle McDonell, Niklas Muennighoff, Chris Ociepa, Jason Phang, Laria Reynolds, Hailey Schoelkopf, Aviya Skowron, Lintang Sutawika, Eric Tang, Anish Thite, Ben Wang, Kevin Wang, and Andy Zou. 2023.  A framework for few-shot language model evaluation.
  * Kim et al. (2022) Young Jin Kim, Rawn Henry, Raffy Fahim, and Hany Hassan Awadalla. 2022.  Who says elephants can’t run: Bringing large scale moe models into cloud scale production.  _arXiv preprint arXiv:2211.10017_.
  * Kwon et al. (2023) Woosuk Kwon, Zhuohan Li, Siyuan Zhuang, Ying Sheng, Lianmin Zheng, Cody Hao Yu, Joseph E. Gonzalez, Hao Zhang, and Ion Stoica. 2023.  Efficient memory management for large language model serving with pagedattention.  In _Proceedings of the ACM SIGOPS 29th Symposium on Operating Systems Principles_.
  * Li et al. (2023) Qingyuan Li, Ran Meng, Yiduo Li, Bo Zhang, Liang Li, Yifan Lu, Xiangxiang Chu, Yerui Sun, and Yuchen Xie. 2023.  A speed odyssey for deployable quantization of llms.  _arXiv preprint arXiv:2311.09550_.
  * Lin et al. (2023) Ji Lin, Jiaming Tang, Haotian Tang, Shang Yang, Xingyu Dang, and Song Han. 2023.  Awq: Activation-aware weight quantization for llm compression and acceleration.  _arXiv preprint arXiv:2306.00978_.
  * Lin et al. (2024) Yujun Lin, Haotian Tang, Shang Yang, Zhekai Zhang, Guangxuan Xiao, Chuang Gan, and Song Han. 2024.  Qserve: W4a8kv4 quantization and system co-design for efficient llm serving.  _arXiv preprint arXiv:2405.04532_.
  * Merity et al. (2016) Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher. 2016.  Pointer sentinel mixture models.  _arXiv preprint arXiv:1609.07843_.
  * Sakaguchi et al. (2021) Keisuke Sakaguchi, Ronan Le Bras, Chandra Bhagavatula, and Yejin Choi. 2021.  Winogrande: An adversarial winograd schema challenge at scale.  _Communications of the ACM_ , 64:99–106.
  * Shao et al. (2023) Wenqi Shao, Mengzhao Chen, Zhaoyang Zhang, Peng Xu, Lirui Zhao, Zhiqian Li, Kaipeng Zhang, Peng Gao, Yu Qiao, and Ping Luo. 2023.  Omniquant: Omnidirectionally calibrated quantization for large language models.  _arXiv preprint arXiv:2308.13137_.
  * Touvron et al. (2023a) Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al. 2023a.  Llama: Open and efficient foundation language models.  _arXiv preprint arXiv:2302.13971_.
  * Touvron et al. (2023b) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al. 2023b.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_.
  * Wei et al. (2023) Xiuying Wei, Yunchen Zhang, Yuhang Li, Xiangguo Zhang, Ruihao Gong, Jinyang Guo, and Xianglong Liu. 2023.  Outlier suppression+: Accurate quantization of large language models by equivalent and optimal shifting and scaling.  _arXiv preprint arXiv:2304.09145_.
  * Xiao et al. (2023) Guangxuan Xiao, Ji Lin, Mickael Seznec, Hao Wu, Julien Demouth, and Song Han. 2023.  Smoothquant: Accurate and efficient post-training quantization for large language models.  In _International Conference on Machine Learning_ , pages 38087–38099. PMLR.
  * Yuan et al. (2023) Zhihang Yuan, Lin Niu, Jiawei Liu, Wenyu Liu, Xinggang Wang, Yuzhang Shang, Guangyu Sun, Qiang Wu, Jiaxiang Wu, and Bingzhe Wu. 2023.  Rptq: Reorder-based post-training quantization for large language models.  _arXiv preprint arXiv:2304.01089_.
  * Zellers et al. (2019) Rowan Zellers, Ari Holtzman, Yonatan Bisk, Ali Farhadi, and Yejin Choi. 2019.  Hellaswag: Can a machine really finish your sentence?  _arXiv preprint arXiv:1905.07830_.
  * Zhang et al. (2022) Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, et al. 2022.  Opt: Open pre-trained transformer language models.  _arXiv preprint arXiv:2205.01068_.
  * Zhao et al. (2023) Yilong Zhao, Chien-Yu Lin, Kan Zhu, Zihao Ye, Lequn Chen, Size Zheng, Luis Ceze, Arvind Krishnamurthy, Tianqi Chen, and Baris Kasikci. 2023.  Atom: Low-bit quantization for efficient and accurate llm serving.  _arXiv preprint arXiv:2310.19102_.