An EM algorithm for absolute
continuous bivariate Pareto distribution

arXiv:1608.02199v4 [stat.CO] 18 Mar 2018

Arabin Kumar Dey∗ , Biplab Paul and Debasis Kundu
Abstract:
Recently [3] used EM algorithm to estimate singular Marshall-Olkin bivariate Pareto distribution. We describe absolute continuous version of
this distribution. We study estimation of the parameters by EM algorithm
both in presence and without presence of location and scale parameters.
Some innovative solutions are provided for different problems arised during implementation of EM algorithm. We also address issues related to
bayesian analysis and proposed methods to compute bayes estimate of the
parameters. A real-life data analysis is also shown for illustrative purpose.
Keywords and phrases: Joint probability density function ; Absolute
Continuous Distribution ; Bivariate pareto distribution ; EM algorithm .

1. Introduction
Recently [3] used EM algorithm to estimate singular Marshall-Olkin bivariate
Pareto distribution. There is no paper on formulation and estimation of absolute
continuous version of this distribution. Few recent works ([11],
[12], [14], [18]) include Marshall-Olkin bivariate distribution. Estimation through
EM algorithm can be seen there. But estimation of absolute continuous version
of bivariate pareto through EM algorithm is not straight forward. The problem becomes more complicated when we consider location and scale parameters
in our problem. Usual gradient descent does not work as bivariate likelihood
is a discontinuous function with respect to location and scale parameters. We
suggest a novel way to handle all related computational problems in most efficient way. The distribution can be used as an alternative model for the data
transformed via peak over threshold method. Bivariate Pareto has wide application in modeling data related to finance, insurance, environmental sciences
and internet network. The dependence structure of absolute continuous version
of bivariate Pareto can also be described by well-known Marshall-Olkin copula
[[19], [15]]. The methodology proposed in this paper works for moderately large
sample.
[5] proposed an absolute continuous bivariate exponential distribution. Later
[13] introduced an absolute continuous bivariate generalized exponential distribution, whose marginals are not generalized exponential distributions. Many
papers related to multivariate extreme value theory talk about absolute continuous/ singular bivariate and multivariate Pareto distribution [[26], [27], [23], [20]].
Development of multivariate Pareto is an asymptotic process in extreme value
∗ Corresponding

Author
1

Dey et al./A sample document

2

theory. None of the paper mentioned specifically anything about MarshalOlkin form of the distribution or its estimation through EM algorithm in
presence of location and scale parameters. A novel EM cum Gradient descent is
the major contribution of this paper. However we address few modifications in
implementation of EM even in case of three parameter set up. The paper contains mix of several issues and implementational issues which together provides
our final proposal algorithm. We show that our proposed algorithm works quite
well for moderately large sample size.
We arrange the paper in the following way. In section 2 we keep the formulation of Marshall-Olkin bivariate Pareto and its absolute continuous version. In
section 3, we describe the EM algorithm for singular case. Different extension
of EM algorithm for absolute continuous version is available in section 4. Some
simulation results show the performance of the algorithm in section 5. In section
6 we show data analysis. Finally we conclude the paper in section 7.
2. Formulation of Marshal-Olkin bivariate pareto
A random variable X is said to have Pareto of second kind, i.e. X ∼ P a(II)(µ, σ, α)
if it has the survival function
F̄X (x; µ, σ, α) = P (X > x) = (1 +

x − µ −α
)
σ

and the probability density function (pdf)
f (x; µ, σ, α) =

α
x − µ −α−1
(1 +
)
σ
σ

with x > µ ∈ R, σ > 0 and α > 0.
Let U0 , U1 and U2 be three independent univariate pareto distributions
P a(II)(0, 1, α0 ), P a(II)(µ1 , σ1 , α1 ) and P a(II)(µ2 , σ2 , α2 ).
We define X1 = min{σ1 U0 + µ1 , U1 } and X2 = min{σ2 U0 + µ2 , U2 }. We can
show that (X1 , X2 ) jointly follow Marshall-Olkin bivariate Pareto distribution
of second kind and we denote it as BV P A(µ1 , σ1 , µ2 , σ2 , α0 , α1 , α2 ).
The joint distribution can be given by

x1 −µ1
x2 −µ2

f1 (x1 , x2 ) if σ1 < σ2
1
2
f (x1 , x2 ) = f2 (x1 , x2 ) if x1σ−µ
> x2σ−µ
1
2


x1 −µ1
x2 −µ2
f0 (x)
if σ1 = σ2 = x
where
f1 (x1 , x2 ) =

α1
(x2 − µ2 ) −(α0 +α2 +1)
(x1 − µ1 ) −(α1 +1)
(α0 + α2 )(1 +
)
(1 +
)
σ1 σ2
σ2
σ1

f2 (x1 , x2 ) =

α2
x1 − µ1 −(α0 +α1 +1)
x2 − µ2 −(α2 +1)
(α0 + α1 )(1 +
)
(1 +
)
σ1 σ2
σ1
σ2

Dey et al./A sample document

3

(a) ξ1

(b) ξ2

(c) ξ3

(d) ξ4

Fig 1. Density plots of BVPAC for different sets of parameters

f0 (x) =

α0
x − µ1 −(α1 +α2 +α0 +1)
(1 +
)
σ1
σ1

We denote absolute continuous part as BVPAC. Surface and contour plots of
the absolute continuous part of the pdf are shown in Figure-1 and Figure-2 respectively. The following four different sets of parameters provide four subfigures
in each figure. ξ1 : µ1 = 0, µ2 = 0, σ1 = 1, σ2 = 0.5, α0 = 1, α1 = 0.3, α2 = 1.4;
ξ2 : µ1 = 1, µ2 = 2, σ1 = 0.4, σ2 = 0.5, α0 = 2, α1 = 1.2, α2 = 1.4; ξ3 : µ1 =
0, µ2 = 0, σ1 = 1.4, σ2 = 0.5, α0 = 1, α1 = 1, α2 = 1.4; ξ4 : µ1 = 0, µ2 = 0, σ1 =
1.4, σ2 = 0.5, α0 = 2, α1 = 0.4, α2 = 0.5.

Dey et al./A sample document

4

(a) ξ1

(b) ξ2

(c) ξ3

(d) ξ4

Fig 2. Contour plots of the density of BVPAC for different sets of parameters

Dey et al./A sample document

5

2.1. Absolute Continuous Extension
We assume that the random vector (Z1 , Z2 ) follows an absolute continuous
bivariate pareto distribution (as described in [12]).
The joint pdf of BVPAC can be written as
fBV P AC (z1 , z2 )

cf1 (z1 , z2 )



= cf (z ; µ , σ , α )f (z ; µ , σ , α + α ) if
PA 1
1
1
1 PA 2
2
2
0
2
=

cf
(z
,
z
)
2 1 2



= cfP A (z1 ; µ1 , σ1 , α0 + α1 )fP A (z2 ; µ2 , σ2 , α2 ) if

(z1 −µ1 )
σ1

<

(z2 −µ2 )
σ2

(z1 −µ1 )
σ1

>

(z2 −µ2 )
σ2

1 +α2
where c is a normalizing constant and c = α0α+α
and fP A (·) denotes the pdf
1 +α2
of univariate distribution.
Likelihood function for the parameters of the BVPAC can be written as

L(µ1 , µ2 , σ1 , σ2 , α0 , α1 , α2 )
= n ln(α0 + α1 + α2 ) − n ln(α1 + α2 ) + n1 ln α1 + n1 ln(α0 + α2 )
X
z2i − µ2
− n1 ln σ1 − n1 ln σ2 − (α0 + α2 + 1)
ln(1 +
)
σ2
i∈I1

− (α1 + 1)

X
i∈I1

+

z1i − µ1
) − n2 σ1 − n2 ln σ2 + n2 ln α2
ln(1 +
σ1

n2 ln(α0 + α1 ) − (α0 + α1 + 1)

X
i∈I2

− (α2 + 1)

X
i∈I2

ln(1 +

z1i − µ1
)
σ1

z2i − µ2
ln(1 +
)
σ2

The three parameter or likelihood function taking µ1 = µ2 = 0 and σ1 =

Dey et al./A sample document

6

σ2 = 1 corresponding to this pdf can be given by
L(α0 , α1 , α2 )

=
+

(n1 + n2 ) ln(α0 + α1 + α2 ) − (n1 + n2 ) ln(α1 + α2 )
X
X
ln fP A (z2i ; 0, 1, α0 + α2 )
ln fP A (z1i ; 0, 1, α1 ) +
i∈I1

i∈I1

+

X

ln fP A (z1i ; 0, 1, α0 + α1 ) +

i∈I2

X

ln fP A (z2i ; 0, 1, α2 )

i∈I2

(n1 + n2 ) ln(α0 + α1 + α2 ) − (n1 + n2 ) ln(α1 + α2 )
X
ln(1 + z1i ) + n1 ln(α0 + α2 )
+ n1 ln α1 − (α1 + 1)

=

i∈I1

− (α0 + α2 + 1)

X

ln(1 + z2i ) + n2 ln(α0 + α1 )

i∈I1

− (α0 + α1 + 1)

X

ln(1 + z1i ) + n2 ln α2

i∈I2

− (α2 + 1)

X

ln(1 + z2i )

i∈I2

Marginal distribution of Z1 and Z2 can be obtained by routine calculation.
Expressions are as follows :
fZ1 (z1 ; µ1 , σ1 , α0 + α1 )
= cfP A (z1 ; µ1 , σ1 , α0 + α1 ) − c

α0
fP A (z1 ; µ1 , σ1 , α0 + α1 + α2 )
(α0 + α1 + α2 )

fZ2 (z2 ; µ2 , σ2 , α0 + α2 )
= cfP A (z2 ; µ2 , σ2 , α0 + α2 ) − c
respectively, where c =

α0
fP A (z2 ; µ2 , σ2 , α0 + α1 + α2 )
(α0 + α1 + α2 )

(α0 +α1 +α2 )
(α1 +α2 ) .

3. EM-algorithm for singular three parameter BVPA distribution
Now we divide our data into three parts :I0 = {(x1i , x2i ) : x1i = x2i } ,I1 = {(x1i , x2i ) : x1i < x2i } , I2 = {(x1i , x2i ) :
x1i > x2i }
We do not know if X1 is U0 or U1 . We also do not know if X2 is U0 or U2 .
So we (
introduce two new random variables (∆1 ; ∆2 ) (
as
0 if X1 = U0
0 if X2 = U0
∆1 =
and
∆2 =
1 if X1 = U1
2 if X2 = U2
The ∆1 and ∆2 are the missing values of the E-M algorithm. To calculate
the E-step we need the conditional distribution of ∆1 and ∆2 .

Dey et al./A sample document

7

Using the definition of X1 , X2 and ∆1 , ∆2 we have :
B For group I0 , both ∆1 and ∆2 are known,
∆1 = ∆2 = 0
B For group I1 , ∆1 is known,∆2 is unknown,
∆1 = 1, ∆2 = 0 or 2
Therefore we need to find out u1 = P (∆2 = 0|I1 )and u2 = P (∆2 = 2|I1 )
B For group I2 , ∆2 is known, ∆1 is unknown,
∆1 = 0 or 1, ∆2 = 2
Moreover we need w1 = P (∆1 = 0|I2 ) and w2 = P (∆1 = 1|I2 )
Since, each posterior probability corresponds to one of the ordering from
Table-3, we calculate u1 , u2 , w1 , w2 using the probability of appropriate ordering.
Ordering
U0 < U1 < U2
U0 < U2 < U1
U1 < U0 < U2
U1 < U2 < U0
U2 < U0 < U1
U2 < U1 < U0

(X1 , X2 )
(U0 , U0 )
(U0 , U0 )
(U1 , U0 )
(U1 , U2 )
(U0 , U2 )
(U1 , U2 )

Group
I0
I0
I1
I1
I2
I2

We have the following expressions for u1 , u2 , w1 , w2
u1 =

P (U1 < U0 < U2 )
P (U1 < U0 < U2 ) + P (U1 < U2 < U0 )

u2 =

P (U1 < U2 < U0 )
P (U1 < U0 < U2 ) + P (U1 < U2 < U0 )

w1 =

P (U2 < U0 < U1 )
P (U2 < U0 < U1 ) + P (U2 < U1 < U0 )

w2 =

P (U2 < U1 < U0 )
P (U2 < U0 < U1 ) + P (U2 < U1 < U0 )

Now we have
Z

∞

[1 − (1 + x)−α1 ]α0 (1 + x)−(α0 +1) (1 + x)−(α2 )

P (U1 < U0 < U2 ) =
0

Z
= α0

∞

(1 + x)−(α0 +α2 +1) − (1 + x)−(α0 +α1 +α2 +1)

0

=

α0 α1
(α0 + α2 + 1)(α0 + α1 + α2 + 1)

Using this above result we evaluate the other probabilites to get values of
0
2
0
1
u1 , u2 , w1 , w2 as : u1 = α0α+α
and u2 = α0α+α
w1 = α0α+α
and w2 = α0α+α
2
2
1
1

Dey et al./A sample document

8

3.1. Pseudo-likelihood expression
We define n0 , n1 , n2 as : n0 = |I0 |, n1 = |I1 |, n2 = |I2 |. where |Ij | for j = 0, 1, 2
denotes the number of elements in the set Ij . Now the pseudo log-likelihood can
be written down as
X
X
X
ln(1 + x2i ))
ln(1 + x1i ) +
ln(1 + xi ) +
L(α0 , α1 , α2 ) = −α0 (
i∈I0

i∈I1

i∈I2

(n0 + u1 n1 + w1 n2 ) ln α0 − α1 (

+

X

ln(1 + xi ) +

(n1 + w2 n2 ) ln α1 − α2 (

X

ln(1 + xi ) +

X

ln(1 + x2i ))

i∈I1 ∪I2

i∈I0

+

ln(1 + x1i ))

i∈I1 ∪I2

i∈I0

+

X

(n2 + u2 n1 ) ln α2

(1)

Therefore M-step involves maximizing (1) with respect to α0 , α1 , α2 at
(t)

(t+1)

α̂0

(t)

n0 + u1 n1 + w1 n2
P
P
.
i∈I0 ln(1 + xi ) +
i∈I2 ln(1 + x1i ) +
i∈I1 ln(1 + x2i )

=P

(2)

(t)

(t+1)

α̂1

=P

(n1 + w2 n2 )
P
ln(1
+
x
i) +
i∈I0
i∈(I1 ∪I2 ) ln(1 + x1i )

(3)

(t)

(t+1)

α̂2

n2 + u2 n1
P
= P
( i∈I0 ln(1 + xi ) + i∈(I1 ∪I2 ) ln(1 + x2i ))

(4)

Therefore the algorithm can be given as
Algorithm 1 EM procedure for bivariate Pareto distribution
1: while ∆Q/Q < tol do
(i)
(i)
(i)
(i)
(i)
(i)
(i)
2:
Compute u1 , u2 , w1 , w2 from α0 , α1 , α2 .
(i+1)

(i+1)

(i+1)

3:
Update α0
, α1
, α2
using Equation (2), (3) and (4).
4:
Calculate Q for the new iterate.
5: end while

4. EM-algorithm for absolute continuous case
We adopt similar existing methods described in [12], we replace the missing
observations n0 falling within I0 and each observation U0 by its estimate n˜0 and
E(U0 |U0 < min{U1 , U2 }). Note that in this case, n0 is a random variable which
1 +α2
has negative binomial with parameters (n1 + n2 ) and α0α+α
.
1 +α2
α0
Clearly n˜0 = (n1 + n2 ) α1 +α2 , E(U0 |U0 < min{U1 , U2 }) = (α0 +α11+α2 −1)
An important restriction for this approximation is that we have to
ensure α0 + α1 + α2 > 1. The restriction will ensure the existence of
the above expectation.

Dey et al./A sample document

9

Therefore under above mentioned restriction, we write the pseudo likelihood
function by replacing the missing observations with its expected value.
X
X
ln(1 + x2i ))
ln(1 + x1i ) +
L(α0 , α1 , α2 ) = −α0 (n˜0 ln(1 + a) +
i∈I1

i∈I2

(n˜0 + u1 n1 + w1 n2 ) ln α0 − α1 (n˜0 ln(1 + a) +

+

X

ln(1 + x1i ))

i∈I1 ∪I2

(n1 + w2 n2 ) ln α1 − α2 (n˜0 ln(1 + a) +

+

X

ln(1 + x2i ))

i∈I1 ∪I2

+

(n2 + u2 n1 ) ln α2

(5)

Therefore M-step involves maximizing (5) with respect to α0 , α1 , α2 at
(t)

(t+1)

α̂0

(t)

n˜ + u1 n1 + w1 n2
P0
P
.
n˜0 ln(1 + a) + i∈I2 ln(1 + x1i ) + i∈I1 ln(1 + x2i )

=

(6)

(t)

(t+1)

α̂1

=

(n1 + w2 n2 )
P
n˜0 ln(1 + a) + i∈(I1 ∪I2 ) ln(1 + x1i )

(7)

(t)

(t+1)

α̂2

=

n2 + u2 n1
P
(n˜0 ln(1 + a) + i∈(I1 ∪I2 ) ln(1 + x2i ))

(8)

Therefore the modified EM algorithm steps would be according to Algorithm2.
Algorithm 2 Modified EM procedure for three parameter absolute continuous
bivariate Pareto distribution
1: while ∆Q/Q < tol do
(i)
(i)
(i)
(i)
(i)
(i)
(i)
2:
Compute u1 , u2 , w1 , w2 , n˜0 , a from α0 , α1 , α2 .
(i+1)

(i+1)

(i+1)

3:
Update α0
, α1
, α2
using Equation (6), (7) and (8).
4:
Calculate Q for the new iterate.
5: end while

Important Remark : The above method fails when α0 + α1 + α2 > 1.
4.1. Proposed Modification :
To make the algorithm valid for any range of parameter, instead of estimating
U0 we estimate ln(1 + U0 ) conditional on U0 < min{U1 , U2 }. Now since ln(x) is
an increasing function in x, we condition on the equivalent event ln(1 + U0 ) <
min{ln(1 + U1 ), ln(1 + U2 )}. Therefore we replace the unknown missing information n0 and ln(1 + U0 ) by
n˜∗ 0 = (n1 + n2 )

α0
α1 + α2

Dey et al./A sample document

10

and
a∗ = E(ln(1 + U0 )| ln(1 + U0 ) < min{ln(1 + U1 ), ln(1 + U2 )}) =

1
(α0 + α1 + α2 )

respectively.
Therefore M-step involves maximizing (5) with respect to α0 , α1 , α2 at
(t+1)

α̂0

=

n˜∗ 0 a∗ +

P

(t)
(t)
n˜∗0 + u1 n1 + w1 n2
P
.
i∈I2 ln(1 + x1i ) +
i∈I1 ln(1 + x2i )

(9)

(t)

(t+1)

α̂1

=

(n1 + w2 n2 )
P
n˜∗ 0 a∗ + i∈(I1 ∪I2 ) ln(1 + x1i )

(10)

(t)

(t+1)

α̂2

=

n2 + u2 n1
P
(n˜∗ 0 a∗ + i∈(I1 ∪I2 ) ln(1 + x2i ))

(11)

Therefore the final EM steps for three parameters can be kept according to
Algorithm-3.
Algorithm 3 Final modified EM procedure for three parameter absolute continuous bivariate Pareto distribution
1: while ∆Q/Q < tol do
(i)
(i)
(i)
(i)
(i)
(i)
(i)
2:
Compute u1 , u2 , w1 , w2 , n˜∗ 0 , a∗ from α0 , α1 , α2 .
(i+1)

(i+1)

(i+1)

3:
Update α0
, α1
, α2
using Equation (9), (10) and (11).
4:
Calculate Q for the new iterate.
5: end while

4.2. Algorithm for seven parameter absolute continuous Pareto
distribution
For seven parameter case we first estimate the location parameters from the
marginal distributions and keep it fixed. Estimates of location parameters can
be replaced by minimum of the marginals of the data. Our EM algorithm update
scale parameters and α0 , α1 , α2 at every iteration. We use one step ahead
gradient descend of σ1 and σ2 based on likelihood of its marginal density. We
can not use the bivariate likelihood function in gradient descend updates of σ1
and σ2 as the likelihood function may not be differentiable with respect to σ1
and σ2 .
Once we get an update of the scale parameters, we first fix I1 and I2 and try
to mimic the EM steps for α0 , α1 and α2 discussed in previous section. We define
z1i = (x1iσˆ−1 µˆ1 ) and z2i = (x2iσˆ−2 µˆ2 ) . If estimates of location and scale parameters
are correct (Z1i , Z2i ) follow a bivariate Pareto distribution with parameter sets
(0, 0, 1, 1, α0 , α1 , α2 ). But one step ahead scale parameter won’t be a correct
estimate of scale parameters, however we assume approximate distribution of
(Z1i , Z2i ) would be Pareto with parameter sets (0, 0, 1, 1, α0 , α1 , α2 ).

Dey et al./A sample document

11

Therefore M-step expression would be same as before :
(t)

(t+1)

α̂0

(t)

n0 + u1 n1 + w1 n2
P
P
.
i∈I0 ln(1 + z1i ) +
i∈I2 ln(1 + z1i ) +
i∈I1 ln(1 + z2i )

=P

(12)

(t)

(n1 + w2 n2 )
P
ln(1
+
z
1i ) +
i∈I0
i∈(I1 ∪I2 ) ln(1 + z1i )

(t+1)

α̂1

=P

(13)

(t)

(t+1)

α̂2

n2 + u2 n1
P
= P
( i∈I0 ln(1 + z1i ) + i∈(I1 ∪I2 ) ln(1 + z2i ))

(14)

We make similar modifications of the above expressions. As usual observations within
I0 and cardinality of I0 would be unknown in this case. We
P
estimate i∈I0 ln(1 + zi ) by n∗0 a∗0 , where a∗0 = E(log(1 + U0 )| log(1 + U0 ) <
0
.
min{log(1 + U1 ), log(1 + U2 )}) = (α0 +α11 +α2 ) and n∗0 by n˜∗0 = (n1 + n2 ) α1α+α
2
Therefore final M-step expression would be same as before :
(t+1)

α̂0

=

n˜∗0 a∗0 +

(t)
(t)
n˜∗0 + u1 n1 + w1 n2
.
P
i∈I2 ln(1 + z1i ) +
i∈I1 ln(1 + z2i )

P

(15)

(t)

(t+1)

α̂1

=

(n1 + w2 n2 )
P
n˜∗0 a∗0 + i∈(I1 ∪I2 ) ln(1 + z1i )

(16)

(t)

(t+1)

α̂2

=

n2 + u2 n1
P
(n˜∗0 a∗0 + i∈(I1 ∪I2 ) ln(1 + z2i ))

(17)

The key intuition of this algorithm is very similar to stochastic gradient
descent. As we increase the number of iteration, EM steps try to force the three
parameters α0 , α1 , α2 to pick up the right direction starting from any value and
gradient descent steps of σ1 and σ2 gradually ensure them to roam around the
actual values. One of the drawback with this approach is that it considers too
many approximations. However this algorithm works even for moderately large
sample sizes. Sometimes it takes lot of time to converge or roam around the
actual value for some really bad sample. Probability of such events are very
low. We stop the calculation after 2000 iteration in such cases. Our numerical
result in next section shows that this does not make any significant effect in the
calculation of mean square error.
Therefore algorithmic steps of our proposed algorithm for seven parameter absolute continuous bivariate Pareto distribution can be kept according to
Algorithm-4.
5. Numerical Results
We use package R 3.2.2 to perform the estimation procedure. All the programs
will be available to author on request. We generate samples of different sizes

Dey et al./A sample document

12

Algorithm 4 EM procedure for seven parameter absolute continuous bivariate
Pareto distribution
1: Take minimum of the marginals as estimates of the location parameters. Initialize other
parameters.
2: while ∆Q/Q < tol do
3:
Compute one step ahead gradient descent update of scale parameters based on likelihood function of the marginals.
4:
Fix I1 and I2 based on the estimated location and scale parameters.
(i)
(i)
(i)
(i)
(i)
(i)
(i)
5:
Compute u1 , u2 , w1 , w2 , n˜∗0 , a∗0 from α0 , α1 , α2 .
(i+1)

(i+1)

(i+1)

using Equation (15), (16) and (17).
, α2
, α1
6:
Update α0
7:
Calculate Q for the new iteration.
8: end while

to calculate the average estimates based on 1000 simulations. We also calculate
mean square error and parametric bootstrap confidence interval for all parameters.
The result shown in Table-1 deals only with three unknown parameters following Algorithm-3. We start EM algorithm taking initial guesses for the parameters as α0 = 2, α1 = 2.2, α2 = 2.4. We try other initial guesses, but average
estimates and MSEs are same. Estimates are calculated based on sample sizes
n = 50, 150, 250, 350, 450 respectively. We use stopping criteria as absolute value
of likelihood changes with respect to previous likelihood at each iteration. In
stopping criteria we use tolerance value as 0.00001. The results shown here are
based on modified algorithm which works for any choice of parameters within
its usual range. The average estimates (AE), mean squared errors (MSE) are
obtained based on 1000 replications.
Table-2 provides the estimates in presence of location and scale parameters
following Algorithm-4. In this case, we take sample size n = 450, 550, 1000, 1500.
However the algorithm works even for smaller sample sizes, although mean
square errors are little higher. Since EM algorithm starts after plug-in the estimates of location parameters as minimum of the maginals, we take initial
values of other parameters as α0 = 2, α1 = 1.2, α2 = 1.4, σ1 = 0.6, σ2 = 0.7.
We find average estimates are closer to true values of the parameters. However
parametric bootstrap confidence interval does not work for location and scale
parameters as it never contains the true parameters.

Dey et al./A sample document

n
50

150

250

350

450

13

parameters

α0

α1

α2

(AE)
MSE
Parametric Bootstrap
Confidence Interval (CI)
(AE)
MSE
Parametric Bootstrap
Confidence Interval (CI)
(AE)
MSE
Parametric Bootstrap
Confidence Interval (CI)
(AE)
MSE
Parametric Bootstrap
Confidence Interval (CI)
(AE)
MSE
Parametric Bootstrap
Confidence Interval (CI)

1.634
0.594
[0.561, 2.733]

0.678
0.278
[0.0003, 1.479]

0.807
0.370
[0.0004, 1.759]

1.865
0.233
[1.023, 2.711]

0.4944abs
0.089
[0.018, 1.066]

0.624
0.125
[0.021, 1.267]

1.9227
0.125
[1.186, 2.630]

0.467
0.054
[0.044, 0.948]

0.577
0.074
[0.047, 1.150]

1.935
0.09
[1.371, 2.595]

0.428
0.033
[0.064, 0.800]

0.543
0.050
[0.087, 0.976]

1.966
0.06
[1.505, 2.486]

0.427
0.027
[0.121, 0.736]

0.528
0.039
[0.138, 0.898]

Table 1
The average estimates (AE), the Mean Square Error (MSE) and parametric bootstrap
confidence interval for α0 = 2, α1 = 0.4 and α2 = 0.5

parameters

µ1

µ2

σ1

AE
MSE

0.101
0.0000038
α0
1.820
0.09738
µ1
0.10117
0.0000026
α0
1.820
0.095
µ1
0.1006574
0.000000848
α0
1.855
0.0655
µ1
0.1004
0.00000039
α0
1.8892
0.038

0.101
0.0000023
α1
0.458
0.02546
µ2
0.10091
0.0000016
α1
0.466
0.0261
µ2
0.1004984
0.00000050
α1
0.481
0.0260
µ2
0.1003
0.000000215
α1
0.4790
0.033

0.753
0.01203
α2
0.626
0.0716
σ1
0.762
0.012
α2
0.632
0.0716
σ1
0.817
0.0085
α2
0.631
0.06122
σ1
0.8095
0.008423629
α2
0.6136
0.0302

n
450
450

550

AE
MSE
parameters
AE
MSE

550
AE
MSE

σ2
0.804
0.0159

σ2
0.807
0.014

1000
AE
MSE
1000

1500

AE
MSE
parameters
AE
MSE

1500
AE
MSE

σ2
0.793
0.012

σ2
0.8213
0.008

Table 2
The average estimates (AE), the Mean Square Error (MSE) for µ1 = 0.1, µ2 = 0.1,
σ1 = 0.8, σ2 = 0.8, α0 = 2, α1 = 0.4 and α2 = 0.5

5.1. Data Set 1 :
The data set is taken from https://archive.ics.uci.edu/ml/machine-learningdatabases. The age of abalone is determined by cutting the shell through the
cone, staining it, and counting the number of rings through a microscope. The
data set contains related measurements. We extract a part of the data for bivariate modeling. We consider only measurements related to female population

Dey et al./A sample document

(a) ξ1

14

(b) ξ2

Fig 3. Survival plots for two marginals of the transformed dataset

Fig 4. Two dimensional density plots for the transformed dataset

where one of the variable is Length as Longest shell measurement and other variable is Diameter which is perpendicular to length. We use peak over threshold
method on this data set.
From [7], we know that peak over threshold method on random variable
U provides polynomial generalized Pareto distribution for any x0 with 1 +
log(G(x0 )) ∈ (0, 1) i.e. P (U > tx0 |U > x0 ) = t−α , t ≥ 1 where G(·) is
the distribution function of U . We choose appropriate t and x0 so that data
should behave more like Pareto distribution. The transformed data set does
not have any singular component. Therefore one possible assumption can be
absolute continuous Marshall Olkin bivariate Pareto.
We also verify our assumption by plotting empirical two dimensional density
plot in Figure-4 which resembles closer to the surface of Marshall-Olkin bivariate
Pareto distribution.
The estimated parameters can be given as µ1 = 10.855, µ2 = 8.632, σ1 =
2.124, σ2 = 1.711, α0 = 3.124, α1 = 1.743 and α2 = 1.602. We can also calculate
the mean square error by parametric bootstrap method. We calculate the MSEs

Dey et al./A sample document

15

based on 1000 bootstrap sample. The corresponding MSE’s are 8.129 ×10−6 ,
5.759 ×10−6 , 0.632, 0.380, 0.759, 0.598, 0.492 respectively.
6. Conclusion:
We observe successful implementation of EM algorithm for absolute continuous bivariate Pareto distribution. The approximation process runs in multiple
stages. Therefore the algorithm works better for larger sample size. Estimation
procedure even in case of three parameters is not straightforward. The paper
also shows some innovative approach to handle the estimation in case of location and scale parameters. But that’s not a problem as it will roam around
closer to the actual value. Since the bivariate likelihood function is discontinuous with respect to location and scale parameters, construction of asymptotic
confidence interval becomes difficult. However parametric bootstrap confidence
interval provides useful solution for this problem. The work of this paper can
be further used for discrimination of several models. The work is in progress.
References
[1] Arnold, B. C. (1967). A note on multivariate distributions with specified
marginals. Journal of the American Statistical Association 62 (320), 1460–
1461.
[2] Asimit, A. V., E. Furman, and R. Vernic (2010). On a multivariate pareto
distribution. Insurance: Mathematics and Economics 46 (2), 308–316.
[3] Asimit, A. V., E. Furman, and R. Vernic (2016). Statistical inference for a
new class of multivariate pareto distributions. Communications in StatisticsSimulation and Computation 45 (2), 456–471.
[4] Balakrishnan, N., R. C. Gupta, D. Kundu, V. Leiva, and A. Sanhueza (2011).
On some mixture models based on the birnbaum–saunders distribution and
associated inference. Journal of Statistical Planning and Inference 141 (7),
2175–2190.
[5] Block, H. W. and A. Basu (1974). A continuous, bivariate exponential extension. Journal of the American Statistical Association 69 (348), 1031–1037.
[6] Dey, A. K. and D. Kundu (2012). Discriminating between bivariate generalized exponential and bivariate weibull distributions. Chilean Journal of
Statistics 3 (1), 93–110.
[7] Falk, M. and A. Guillou (2008). Peaks-over-threshold stability of multivariate generalized pareto distributions. Journal of Multivariate Analysis 99 (4),
715–734.
[8] Gupta, R. D. and D. Kundu (1999). Theory & methods: Generalized exponential distributions. Australian & New Zealand Journal of Statistics 41 (2),
173–188.
[9] Khosravi, M., D. Kundu, and A. Jamalizadeh (2015). On bivariate and a
mixture of bivariate birnbaum–saunders distributions. Statistical Methodology 23, 1–17.

Dey et al./A sample document

16

[10] Kundu, D. (2012). On sarhan-balakrishnan bivariate distribution. Journal
of Statistics Applications & Probability 1 (3), 163.
[11] Kundu, D. and R. D. Gupta (2009). Bivariate generalized exponential
distribution. Journal of Multivariate Analysis 100 (4), 581–593.
[12] Kundu, D. and R. D. Gupta (2010). A class of absolutely continuous bivariate distributions. Statistical Methodology 7 (4), 464–477.
[13] Kundu, D. and R. D. Gupta (2011). Absolute continuous bivariate generalized exponential distribution. AStA Advances in Statistical Analysis 95 (2),
169–185.
[14] Kundu, D., A. Kumar, and A. K. Gupta (2015). Absolute continuous
multivariate generalized exponential distribution. Sankhya B 77 (2), 175–206.
[15] Marshall, A. W. (1996). Copulas, marginals, and joint distributions. Lecture
Notes-Monograph Series, 213–222.
[16] Marshall, A. W. and I. Olkin (1967). A multivariate exponential distribution. Journal of the American Statistical Association 62 (317), 30–44.
[17] McLachlan, G. and D. Peel (2004). Finite mixture models. John Wiley &
Sons.
[18] Mirhosseini, S. M., M. Amini, D. Kundu, and A. Dolati (2015). On a new
absolutely continuous bivariate generalized exponential distribution. Statistical Methods & Applications 24 (1), 61–83.
[19] Nelsen, R. B. (2007). An introduction to copulas. Springer Science &
Business Media.
[20] Rakonczai, P. and A. Zempléni (2012). Bivariate generalized pareto distribution in practice: models and estimation. Environmetrics 23 (3), 219–227.
[21] Ristić, M. M. and D. Kundu (2015). Marshall-olkin generalized exponential
distribution. Metron 73 (3), 317–333.
[22] Sarhan, A. M. and N. Balakrishnan (2007). A new class of bivariate distributions and its mixture. Journal of Multivariate Analysis 98 (7), 1508–1527.
[23] Smith, R. L. (1994). Multivariate threshold methods. In Extreme Value
Theory and Applications, pp. 225–248. Springer.
[24] Smith, R. L., J. A. Tawn, and S. G. Coles (1997). Markov chain models
for threshold exceedances. Biometrika 84 (2), 249–268.
[25] Tsay, R. S. (2014). An introduction to analysis of financial data with R.
John Wiley & Sons.
[26] Yeh, H.-C. (2000). Two multivariate pareto distributions and their related
inferences. Bulletin of the Institute of Mathematics, Academia Sinica. 28 (2),
71–86.
[27] Yeh, H.-C. (2004). Some properties and characterizations for generalized
multivariate pareto distributions. Journal of Multivariate Analysis 88 (1),
47–60.
Department of Mathematics,
IIT Guwahati
e-mail: arabin.k.dey@gmail.com

Department of Mathematics,
IIT Guwahati
e-mail: biplab.paul@iitg.ac.in

Dey et al./A sample document
Department of Mathematics and Statistics
IIT Kanpur
Kanpur, India
e-mail: kundu@iitk.ac.in

17

