Ibrahim Elbatal
email: aiielbatal@imamu.edu.sa

Farrukh Jamal
email: bfarrukhjamalmphil@gmail.com

Christophe Chesneau
email: christophe.chesneau@unicaen.fr

Mohammed Elgarhy
email: dm_elgarhy85@yahoo.com

Sharifah Alrajhi
email: esaalrajhi@kau.edu.sa

The Modified Beta Gompertz Distribution: Theory and Applications

Keywords: Modified beta generator, Gompertz distribution, Maximum likelihood estimation. 2000 MSC: 60E05, 62E15, 62F10

In this paper, we introduce a new continuous probability distribution with five parameters called the modified beta Gompertz distribution. It is derived from the modified beta generator proposed by Nadarajah, Teimouri and Shih (2014) and the Gompertz distribution. By investigating its mathematical and practical aspects, we prove that it is quite flexible and can be used effectively in modeling a wide variety of real phenomenon. Among others, we provide useful expansions of crucial functions, quantile function, moments, incomplete moments, moment generating function, entropies and order statistics. We discuss the estimation of the model parameters by the obtained maximum likelihood method. We also present a simulation study in order to test validity of maximum likelihood estimators. Finally, we illustrate the flexibility of the distribution by means of two applications to real data.

Introduction

The Gompertz distribution was initially introduced by [START_REF] Gompertz | On the nature of the function expressive of the law of human mortality and on the new mode of determining the value of life contingencies[END_REF] to describe human mortality and provide actuarial tables. The literature about the use of the Gompertz distribution in applied areas is enormous. A nice review can be found in [START_REF] Tjørve | The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family[END_REF], and the references there in. From a mathematical point of view, the cumulative probability density function (cdf) of the Gompertz distribution with parameters λ > 0 and α > 0 is given by G(x) = 1 -e -λ α (e αx -1) , x > 0.

The related probability density function (pdf) is given by g(x) = λe αx e -λ α (e αx -1) , x > 0.

It can be viewed as a generalization of the exponential distribution (obtained with α → 0) and thus an alternative to the gamma or Weibull distribution. A feature of the Gompertz distribution is that g(x) is unimodal and has positive skewness whereas the related hazard rate function (hrf) given by h(x) = g(x)/(1 -G(x)) is increasing. In order to increase the flexibility of the Gompertz distribution, further extensions have been proposed. A natural one is the generalized Gompertz distribution introduced by El- [START_REF] El-Gohary | The generalized Gompertz distribution[END_REF]. By introducing an exponent parameter a > 0, the related cdf is given by

F (x) = 1 -e -λ
α (e αx -1) a , x > 0.

The related applications show that a plays an important role in term of model flexibility. This idea was then extended by [START_REF] Jafari | The beta-Gompertz distribution[END_REF] by using the so-called beta generator introduced by [START_REF] Eugene | Beta-normal distribution and its applications[END_REF]. The related cdf is given by

F (x) = 1 B(a, b)
1-e -λ α (e αx -1) 0 t a-1 (1 -t) b-1 dt = I

1-e -λ α (e αx -1) (a, b), x > 0.

(1)

where a, b > 0, B(a, b) denotes the beta function defined by B(a, b) = 1 0 t a-1 (1t) b-1 dt and I x (a, b) denotes the incomplete beta function ratio defined by I x (a, b) = (1/B(a, b))

x 0 t a-1 (1 -t) b-1 dt, x ∈ [0, 1]. This distribution has been recently extended by Benkhelifa (2017) with a five parameters distribution. It is based on the beta generator and the generalized Gompertz distribution.

Motivated by the emergence of complex data from many applied areas, other extended Gompertz distributions have been proposed in the literature. See for instance, [START_REF] El-Damcese | Generalized Exponential Gompertz Distribution[END_REF] who consider the Odd Generalized Exponential generator introduced by [START_REF] Tahir | The Odd Generalized Exponential Family of Distributions with Applications[END_REF], [START_REF] Roozegar | The McDonald Gompertz distribution: properties and applications[END_REF] who use the the McDonald generator introduced by [START_REF] Alexander | Generalized beta-generated distributions[END_REF], [START_REF] Moniem | Transmuted Gompertz Distribution[END_REF] and [START_REF] Khan | Transmuted Generalized Gompertz distribution with application[END_REF] who apply the transmuted generator introduced by Shaw and [START_REF] Buckley | The alchemy of probability distributions: beyond Gram-Charlier expansions and a skewkurtotic-normal distribution from a rank transmutation map[END_REF], [START_REF] Lima | The Kumaraswamy Gompertz distribution[END_REF] and [START_REF] Chukwu | On kumaraswamy gompertz makeham distribution[END_REF] who use the Kumaraswamy generator, and Benkhelifa (2017), [START_REF] Yaghoobzadeh | A new generalization of the Marshall-Olkin Gompertz distribution[END_REF] who consider the Marshall-Olkin generator introduced by [START_REF] Marshall | A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families[END_REF] and [START_REF] Shadrokh | The Beta Gompertz Geometric Distribution: Mathematical Properties and Applications[END_REF] who consider the Beta-G and Geometric generators.

In this paper, we present and study a distribution with five parameters extending the Gompertz distribution. It is based on the modified beta generator developed by [START_REF] Nadarajah | Modified Beta Distributions[END_REF] (which can also be viewed as a modification of the beta Marshall-Olkin generator developed by [START_REF] Alizadeh | The beta Marshall-Olkin family of distributions[END_REF]). The advantage of this generator is to nicely combine the advantages of the beta generator of [START_REF] Eugene | Beta-normal distribution and its applications[END_REF] and the Marshall-Olkin generator of [START_REF] Marshall | A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families[END_REF]. To the best of our knowledge, its application to the Gompertz distribution has never been considered before. We provide a comprehensive description of its general mathematical properties (expansions of the cdf and pdf, quantile function, various kinds of moments, moment generating function, entropies and order statistics). The estimation of the model parameters by maximum likelihood is then discussed. Finally, we explore applications to real data sets that illustrate the usefulness of the proposed model.

The rest of the paper is organized as follows. Section 2 describes the considered modified beta Gompertz distribution. Some mathematical properties are investigated in Section 3. Section 4 provides the necessary to the estimation of the unknown parameters with the maximum likelihood method. A simulation study is performed in order to test validity of the obtained maximum likelihood estimators. To illustrate the flexibility of the resulting model, applications to two real life data sets are also given.

The Modified Beta Gompertz Distribution

Let c > 0, G(x) be a cdf and g(x) be a related pdf. The modified beta generator introduced by [START_REF] Nadarajah | Modified Beta Distributions[END_REF] is characterized by the cdf given by

F (x) = I cG(x) 1-(1-c)G(x) (a, b), (2)
By differentiation of F (x), a pdf is given by

f (x) = c a g(x) [G(x)] a-1 [1 -G(x)] b-1 B(a, b) [1 -(1 -c)G(x)] a+b , x ∈ R. ( 3
)
The hrf is given by

h(x) = c a g(x) [G(x)] a-1 [1 -G(x)] b-1 B(a, b) [1 -(1 -c)G(x)] a+b 1 -I cG(x) 1-(1-c)G(x) (a, b) , x ∈ R.
Let us now present our main distribution of interest. Using the cdf G(x) of the Gompertz distribution with parameters λ > 0 and α > 0 as baseline, the cdf given by (2) becomes

F (x) = I c 1-e -λ α (e αx -1) 1-(1-c) 1-e -λ α (e αx -1) (a, b), x > 0. ( 4
)
The related distribution will be call the modified beta Gompertz distribution (MBGz distribution for short), also denoted by MBGz(λ, α, a, b, c). The related pdf (3) is given by

f (x) = c a λe αx e -λb α (e αx -1) 1 -e -λ α (e αx -1) a-1 B(a, b) 1 -(1 -c) 1 -e -λ α (e αx -1) a+b , x > 0. ( 5
)
The hrf is given by

h(x) = c a λe αx e -λb α (e αx -1) 1 -e -λ α (e αx -1) a-1 B(a, b) 1 -(1 -c) 1 -e -λ α (e αx -1) a+b       1 -I c 1-e -λ α (e αx -1) 1-(1-c) 1-e -λ α (e αx -1) (a, b)       , x > 0. ( 6
)
Figure 1 shows the plots for f (x) and h(x) for selected parameter values λ, α, a, b, c. We observe that these functions can take various curvature forms depending on the parameter values, showing the increasing of the flexibility of the former Gompertz distribution.

λ = 0.9 α = 1 a = 0.2 b = 1 c = 0.1 λ = 1 α = 0.5 a = 3 b = 0.5 c = 1 λ = 0.1 α = 0.2 a = 0.1 b = 0.2 c = 2 λ = 1 α = 1 a = 0.2 b = 0.2 c = 2
Figure 1: Some plots of the pdf f (x) (a) and some plots for the hrf h(x) (b).

A strong point of the MBGz distribution is to contain different useful distributions in the literature. The most popular of them are listed below.

• When c = 1/(1 -θ) with θ ∈ (0, 1) (θ is a proportion parameter), we obtain the beta Gompertz geometric distribution introduced by Shadrokh and Yaghoobzadeh (2018), i.e. with cdf

F (x) = I 1-e -λ α (e αx -1) 1-θe -λ α (e αx -1) (a, b), x > 0.
However, this distribution excludes the case c ∈ (0, 1) by construction. The importance of small values for c can also be determinant in the applications (see Section 4).

• When c = 1, we get the beta Gompertz distribution with four parameters introduced by [START_REF] Jafari | The beta-Gompertz distribution[END_REF], i.e. with cdf

F (x) = I 1-e -λ
α (e αx -1) (a, b), x > 0.

• When c = b = 1, we get the generalized Gompertz distribution studied by [START_REF] El-Gohary | The generalized Gompertz distribution[END_REF], i.e. with cdf

F (x) = 1 -e -λ
α (e αx -1) a , x > 0.

• When a = b = 1 and c = 1 θ with θ > 1 we get the a particular case of the Marshall-Olkin extended generalized Gompertz distribution introduced by Benkhelifa (2017), i.e. with cdf

F (x) = 1 -e -λ α (e αx -1) θ + (1 -θ) 1 -e -λ
α (e αx -1)

, x > 0.

• When a = b = c = 1, we get the Gompertz distribution introduced by Gompertz (1825), i.e. with cdf

F (x) = 1 -e -λ
α (e αx -1) , x > 0.

• When c = 1 and α → 0, we get beta Exponential distribution studied by Nadarajah and Kotz (2006), i.e. with cdf

F (x) = I 1-e -λx (a, b), x > 0.
• When b = c = 1 and α → 0, we get the generalized exponential distribution studied by [START_REF] Gupta | Generalized Exponential Distributions[END_REF], i.e. with cdf

F (x) = 1 -e -λx a , x > 0.
• When a = b = c = 1 and α → 0 we get the exponential distribution, i.e. with cdf

F (x) = 1 -e -λx , x > 0.
3. Some mathematical properties

On the shapes of the pdf

The shapes of f (x) given by (5) can be described analytically. As usual, the critical points x * of the pdf f (x) satisfies ∂ ∂x ln(f (x * )) = 0, with

∂ ∂x ln(f (x)) = α -bλe αx + (a -1) λe αx e -λ α (e αx -1)
1 -e -λ α (e αx -1)

+ (a + b)(1 -c) λe αx e -λ α (e αx -1) c + (1 -c)e -λ α (e αx -1) . A point x * corresponds to a local maximum if ∂ 2 ∂x 2 ln(f (x * )) < 0, a local minimum if ∂ 2 ∂x 2 ln(f (x * )) > 0 and a point of inflection if ∂ 2 ∂x 2 ln(f (x * )) = 0.
Let us now study the asymptotic properties of f (x). We have

f (x) ∼ c a B(a, b) λ a x a-1 , x → 0.
So, for a ∈ (0, 1), we have lim x→0 f (x) = +∞, for a = 1, we have lim x→0 f (x) = bcλ and for a > 1, we have lim x→0 f (x) = 0. We have

f (x) ∼ 1 c b B(a, b)
λe αx e λb α e -λb α e αx , x → +∞.

Thus lim x→+∞ f (x) = 0 in all cases. Figure 1 (a) illustrates these points for selected parameters.

On the shapes of the hrf

Similarly to the pdf, the critical points x * of the hrf h(x) given by ( 6) satisfies

∂ ∂x ln(h(x * )) = 0, with ∂ ∂x ln(h(x)) = α -bλe αx + (a -1) λe αx e -λ α (e αx -1)
1 -e -λ α (e αx -1)

+ (a + b)(1 -c) λe αx e -λ α (e αx -1) c + (1 -c)e -λ α (e αx -1) + c a λe αx e -λb α (e αx -1) 1 -e -λ α (e αx -1) a-1 B(a, b) 1 -(1 -c) 1 -e -λ α (e αx -1) a+b       1 -I c 1-e -λ α (e αx -1) 1-(1-c) 1-e -λ α (e αx -1) (a, b)       . A point x * corresponds to a local maximum if ∂ 2 ∂x 2 ln(h(x * )) < 0, a local minimum if ∂ 2 ∂x 2 ln(h(x * )) > 0 and a point of inflection if ∂ 2 ∂x 2 ln(h(x * )) = 0. We also have h(x) ∼ c a B(a, b) λ a x a-1 , x → 0.
So, for a ∈ (0, 1), we have lim x→0 h(x) = +∞, for a = 1, we have lim x→0 h(x) = bcλ and for a > 1, we have lim x→0 h(x) = 0. We have

h(x) ∼ bλe αx , x → +∞.
Thus lim x→+∞ h(x) = +∞ in all cases. Figure 1 (b) illustrates these points for selected parameters.

Linear representation

Let us determine useful linear representations for F (x) given by (4) and f (x) given by ( 5). First of all, let us suppose that c ∈ (0, 1). It follows from the generalized binomial formula, i.e. (1

+ z) γ = +∞ k=0 γ k z k for |z| < 1 and γ ∈ R, with γ k = γ(γ -1) . . . (γ -k + 1) k! , that F (x) = 1 B(a, b) cG(x) 1-(1-c)G(x) 0 t a-1 (1 -t) b-1 dt = 1 B(a, b) +∞ k=0 b -1 k (-1) k cG(x) 1-(1-c)G(x) 0 t a+k-1 dt = 1 B(a, b) +∞ k=0 b -1 k (-1) k a + k cG(x) 1 -(1 -c)G(x) a+k .
On the other hand, using again the generalized binomial formula, we obtain

cG(x) 1 -(1 -c)G(x) a+k = c a+k +∞ =0 -(a + k) (-1) (1 -c) [G(x)] +a+k . (7)
In a similar manner, we have

[G(x)] +a+k = 1 -e -λ
α (e αx -1)

+a+k = +∞ m=0 + a + k m (-1) m (1 -H m (x)),
where H m (x) = 1 -e -mλ α (e αx -1) is the cdf of a Gompertz distribution with parameters mλ and α. Combining these equalities, we obtain the following series expansion:

F (x) = +∞ m=0 v m (1 -H m (x)), (8)
where

v m = (-1) m B(a, b) +∞ k=0 +∞ =0 + a + k m -(a + k) b -1 k c a+k (-1) k+ (1 -c) 1 a + k .
By derivation of F (x), f (x) can be expressed as

f (x) = +∞ m=0 w m h m (x), (9)
where w m = -v m and h m (x) is the pdf of a Gompertz distribution with parameters mλ and α.

For the case c > 1, we must do some transformation for the equation ( 7) in order to apply the generalized binomial formula. We can write

cG(x) 1 -(1 -c)G(x) a+k = G(x) 1 -(1 -1 c )(1 -G(x)) a+k = [G(x)] a+k +∞ =0 -(a + k) (-1) (c -1) c -[1 -G(x)] = +∞ =0 q=0 -(a + k) q (-1) +q (c -1) c -[G(x)] q+a+k .
On the other hand, we have

[G(x)] q+a+k = +∞ m=0 q + a + k m (-1) m (1 -H m (x)).
Therefore, we can write F (x) as ( 8) with

v * m = (-1) m B(a, b) +∞ k=0 +∞ =0 q=0 -(a + k) q q + a + k m b -1 k (-1) +q+k (c -1) c -1 a + k ,
and f (x) as ( 9) with w m = -v * m (and still h m (x) is the pdf of a Gompertz distribution with parameters mλ and α). For the sake of simplicity, we shall refer to the form (9) far all series representation of f (x), whatever c ∈ (0, 1) or c > 1.

Hereafter, we denote by X a random variable having the cdf F (x) given by (4) (and the pdf f (x) given by ( 5)) and by Y m a random variable following the Gompertz distribution with parameters mλ and α, i.e. having the cdf H m (x) (and the pdf h m (x)).

Quantile function

The quantile function of X is given by

Q(u) = 1 α ln 1 - α λ ln 1 - I -1 u (a, b) c + (1 -c)I -1 u (a, b)
, u ∈ (0, 1),

where

I -1 u (a, b) denotes the inverse of I u (a, b). It satisfies F (Q(u)) = Q(F (u)) = u. Using Nadarajah et al. (2014), one can show that Q(u) ∼ 1 λc a 1 a B(a, b) 1 a u 1 a , u → 0.
From Q(u), we can simulate the MBGz distribution. Indeed, let U be a random variable following the uniform distribution over (0, 1). Then the random variable X = Q(U ) follows the MBGz distribution.

The median of X is given by M = Q(1/2). We can also use Q(u) to define skewness measures. Let us just introduce the Bowley skewness based on quartiles and the Moors kurtosis respectively defined by

B = Q(3/4) + Q(1/4) -2Q(1/2) Q(3/4) -Q(1/4) , M o = Q(7/8) -Q(5/8) + Q(3/8) -Q(1/8) Q(6/8) -Q(2/8
) .

Contrary to γ 1 and γ 2 , these quantities have the advantage to be always defined. We refer to [START_REF] Kenney | Mathematics of Statistics[END_REF] and [START_REF] Moors | A quantile alternative for kurtosis[END_REF].

Moments

Let r be a positive integer. The r-th ordinary moment of X is defined by µ

r = E (X r ) = +∞ -∞ x r f (x)dx.
Using the linear representation given by ( 9), we can express µ r as

µ r = +∞ m=0 w m +∞ -∞ x r h m (x)dx = +∞ m=0 w m E(Y r m ).
By doing the change of variables u = e αx , we obtain

E(Y r m ) = mλ α r+1 e mλ α +∞ 1 (ln u) r e -mλ α u du.
This integral has connections with the so-called generalized integro-exponential function. Further developments can be found in [START_REF] Milgram | The generalized integro-exponential function[END_REF] and [START_REF] Lenart | The moments of the Gompertz distribution and maximum likelihood estimation of its parameters[END_REF]. Therefore we have

µ r = +∞ m=0 w m mλ α r+1 e mλ α +∞ 1 (ln u) r e -mλ α u du.
Obviously, the mean of X is given by E(X) = µ 1 and the variance of X is given by V(X) = µ 2 -(µ 1 ) 2 .

Skewness

The r-th central moment of X is given by µ r = E [(X -µ 1 ) r ]. It follows from the binomial formula that

µ r = r k=0 r k (-1) k (µ 1 ) k µ r-k .
On the other side, the r-th cumulants of X can be obtained via the equation:

κ r = µ r - r-1 k=1 r -1 k -1 κ k µ r-k , with κ 1 = µ 1 . The skewness of X is given by γ 1 = κ 3 /κ 3/2 2
and the kurtosis of X is given by γ 2 = κ 4 /κ 2 2 . One can also introduce the MacGillivray skewness given by

ρ(u) = Q(1 -u) + Q(u) -2Q(1/2) Q(1 -u) -Q(u) , u ∈ (0, 1).
It illustrates the effects of the parameters a, b, α and λ on the skewness. Further details can be found in MacGillivray (1986).

Moment generating function

The moment generating function of X is given by M X (t) = E e tX = +∞ -∞ e tx f (x)dx. Using (9), we have

M X (t) = +∞ m=0 w m +∞ -∞ e tx h m (x)dx = +∞ m=0 w m M Ym (t),
where M Ym (t) = E(e tYm ), the moment generating function of Y m . Doing successively the change of variables u = e αx and the change of variable v = mλ α u, we obtain

M Ym (t) = mλ α e mλ α +∞ 1 u t α e -mλ α u du = e mλ α α mλ t α +∞ mλ α v t α e -v dv = e mλ α α mλ t α Γ t α + 1, mλ α ,
where Γ(d, x) denotes the complementary incomplete gamma function defined by Γ(d, x) = +∞ x t d-1 e -t dt. Therefore we can write

M X (t) = +∞ m=0 w m e mλ α α mλ t α Γ t α + 1, mλ α .
Alternatively, using the moments of X, one can write

M X (t) = +∞ r=0 t r r! µ r = +∞ r=0 +∞ m=0 t r r! w m mλ α r+1 e mλ α
+∞ 1 (ln u) r e -mλ α u du.

Incomplete moments and mean deviations

The r-th incomplete moment of X is defined by m r (t) = E X r 1 {X≤t} = t -∞ x r f (x)dx. Using (9), we can express m r (t) as

m r (t) = +∞ m=0 w m t -∞ x r h m (x)dx.
Doing successively the change of variables u = e αx , we obtain

t -∞ x r h m (x)dx = mλ α r+1 e mλ α
e αt 1 (ln u) r e -mλ α u du.

The mean deviation of X about the mean is given by

δ 1 = E(|X -µ 1 |) = 2µ 1 F (µ 1 ) -2m 1 (µ 1 ),
where m 1 (t) denote the first incomplete moment. The mean deviation of X about the median M = Q(1/2) is given by

δ 2 = E(|X -M |) = µ 1 -2m 1 (M ).

Entropies

Let us now investigate different kinds of entropy measures. The Rényi entropy of X is defined by a+b) .

I γ (X) = 1 1 -γ ln +∞ -∞ [f (x)] γ dx , with γ > 0 and γ = 1. It follows from (3) that [f (x)] γ = c aγ [g(x)] γ [G(x)] γ(a-1) [1 -G(x)] γ(b-1) B(a, b) γ [1 -(1 -c)G(x)] γ(
The generalized binomial formula implies that -1) .

[G(x)] γ(a-1) [1 -(1 -c)G(x)] γ(a+b) = +∞ k=0 -γ(a + b) k (-1) k (1 -c) k [G(x)] k+γ(a
Similarly, we have

[G(x)] k+γ(a-1) = +∞ =0 k + γ(a -1) (-1) [1 -G(x)] . Therefore [f (x)] γ = c aγ B(a, b) γ +∞ k=0 +∞ =0 -γ(a + b) k k + γ(a -1) (-1) k+ (1 -c) k [1 -G(x)] +γ(b-1) [g(x)] γ .
By doing the change of variable u = e αx and the change of variable v = ( + γb) λ α u, we get

+∞ -∞ [1 -G(x)] +γ(b-1) [g(x)] γ dx = +∞ 0 e -( +γb) λ α (e αx -1) λ γ e αγx dx = λ γ 1 α e ( +γb) λ α +∞ 1 u γ-1 e -( +γb) λ α u du = α γ-1 ( + γb) γ e ( +γb) λ α +∞ ( +γb) λ α v γ-1 e -v dv = α γ-1 ( + γb) γ e ( +γb) λ α Γ γ, ( + γb) λ α .
By putting the above equalities together, we have

I γ (X) = 1 1 -γ αγ ln(c) -γ ln(B(a, b)) + (γ -1) ln(α) + γbλ α + ln +∞ k=0 +∞ =0 -γ(a + b) k k + γ(a -1) (-1) k+ (1 -c) k e λ α ( + γb) γ Γ γ, ( + γb) λ α .
The Shannon entropy of X is defined by S(X) = E(-ln[f (X)]) is a particular case of the Rényi entropy when γ tends to 1 + . The γ-entropy is defined by

H γ (X) = 1 γ -1 ln 1 - +∞ -∞ [f (x)] γ dx .
Using the expansion above, we obtain

H γ (X) = 1 γ -1 ln 1 - c aγ α γ-1 e γb λ α B(a, b) γ × +∞ k=0 +∞ =0 -γ(a + b) k k + γ(a -1) (-1) k+ (1 -c) k e λ α ( + γb) γ Γ γ, ( + γb) λ α .

Order statistics

Let X 1 , . . . , X n be the random sample from X and X i:n be the i-th order statistic. Then the pdf of X i:n is given by

f i:n (x) = n! (i -1)!(n -i)! f (x)[F (x)] i-1 [1 -F (x)] n-i = n! (i -1)!(n -i)! n-i j=0 n -i j (-1) j f (x)[F (x)] j+i-1 .
It follows from ( 8) and ( 9) that

f i:n (x) = n! (i -1)!(n -i)! n-i j=0 n -i j (-1) j +∞ m=0 w m h m (x) +∞ k=0 v k (1 -H k (x)) j+i-1
.

Using a result from [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF], power series raised to a positive power as follows

+∞ k=0 a k x k n = +∞ k=0 d n,k x k ,
where the coefficients (d n,k ) k∈N are determined from the recurrence equation:

d n,0 = a n 0 and, for any m ≥ 1, d n,m = (1/(ma 0 )) m k=1 (k(n + 1) -m)a k d n,m-k . Therefore, noticing that 1 -H k (x) = e -λ α (e αx -1) k , we have +∞ k=0 v k (1 -H k (x)) j+i-1 = +∞ k=0 d j+i-1,k (1 -H k (x)), 0 and, for any m ≥ 1, d j+i-1,m = 1 mv 0 m k=1 (k(j+i)-m)v k d j+i-1,m-k .
By combining the equalities above, we obtain

f i:n (x) = n! (i -1)!(n -i)! n-i j=0 n -i j (-1) j +∞ m=0 +∞ k=0 w m d j+i-1,k h m (x)(1 -H k (x)).
Finally, one can observe that h m (x)(1 -H k (x)) = mλe αx e -(m+k)λ α (e αx -1) = m m+k u m+k (x), where u m+k (x) denotes the pdf of the Gompertz distribution with parameters (m + k)λ and α. So the pdf of i-th order statistic of the MBGz distribution can be expressed as a linear combination of Gompertz pdfs, i.e.

f i:n (x) = n! (i -1)!(n -i)! n-i j=0 n -i j (-1) j +∞ m=0 +∞ k=0 w m d j+i-1,k m m + k u m+k (x).
Let r be a positive integer. Then the r-th ordinary moment of X i:n can be expressed as

E(X r i:n ) = +∞ -∞ x r f i:n (x)dx = n! (i -1)!(n -i)! n-i j=0 n -i j (-1) j +∞ m=0 +∞ k=0 w m d j+i-1,k mλ α r+1 e (m+k)λ α +∞ 1 (ln u) r e -(m+k)λ α u du.

Statistical inference

Maximum likelihood estimation

We now investigate the estimation of the parameters of the MBGz distribution. Let x 1 , . . . , x n be n observed values from the MBGz distribution and ξ = (λ, α, a, b, c) be the vector of unknown parameters. The log likelihood function is given by

(ξ) = an ln(c) + n ln(λ) + α n i=1 x i - λb α n i=1 (e αx i -1) + (a -1) n i=1 ln 1 -e -λ α (e αx i -1) -n ln(B(a, b)) -(a + b) n i=1 ln 1 -(1 -c) 1 -e -λ α (e αx i -1)
.

The maximum likelihood estimators of the parameters are obtained by maximizing the log likelihood function. They can be obtained by solving the non-linear equations:

∂ ∂λ (ξ) = 0, ∂ ∂α (ξ) = 0, ∂ ∂a (ξ) = 0, ∂ ∂b (ξ) = 0, ∂ ∂c (ξ) = 0 with ∂ (ξ) ∂λ = n λ - b α n i=1 (e αx i -1) + (a -1) n i=1 1 α (e αx i -1)e -λ α (e αx i -1)
1 -e -λ α (e αx i -1)

+ (a + b) n i=1 (1 -c) 1 α (e αx i -1)e -λ α (e αx i -1) 1 -(1 -c) 1 -e -λ α (e αx i -1) , ∂ (ξ) ∂α = n i=1 x i - λb α n i=1
x i e αx i -1 α (e αx i -1)

+ (a -1) n i=1 λ α e -λ α (e αx i -1) x i e αx i -1 α (e αx i -1) 1 -e -λ α (e αx i -1) + (a + b) n i=1 (1 -c)e -λ α (e αx i -1) x i e αx i -1 α (e αx i -1) 1 -(1 -c) 1 -e -λ α (e αx i -1) , by setting B (1,0) (a, b) = ∂ ∂a B(a, b) and B (0,1) (a, b) = ∂ ∂b B(a, b) (one can remark that B (1,0) (a, b) = ψ(a) -ψ(a + b) and B (0,1) (a, b) = ψ(b) -ψ(a + b)
, where ψ(x) denotes the so called digamma function),

∂ (ξ) ∂a = n ln c + n i=1 ln 1 -e -λ α (e αx i -1) -n B (1,0) (a, b) B(a, b) - n i=1 ln 1 -(1 -c) 1 -e -λ α (e αx i -1) , ∂ (ξ) ∂b = - λ α n i=1 (e αx i -1) -n B (0,1) (a, b) B(a, b) - n i=1 ln 1 -(1 -c) 1 -e -λ α (e αx i -1)
and

∂ (ξ) ∂c = an c -(a + b) n i=1 1 -e -λ
α (e αx i -1)

1 + (1 -c) 1 -e -λ
α (e αx i -1)

.

We can solve the above non-linear equations simultaneously. Mathematical package can be used to get the maximum likelihood estimators of the unknown parameters. Also, all the second order derivatives exist. As usual, the normal approximation for the maximum likelihood estimators can be used for constructing approximate confidence intervals, confidence regions and testing hypotheses of λ, α, a, b, c.

Simulation

It is very difficult to compare the theoretical performances of the different maximum likelihood estimates (MLEs) for the MBGz distribution. Therefore, simulation is needed to compare the performances of the MLE mainly with respect to their mean square errors for different sample sizes. A numerical study is performed using Mathematica 9 software. Different sample sizes are considered through the experiments at size n = 50, 100 and 150. The experiment will be repeated 3000 times. In each experiment, the estimates of the parameters will be obtained by maximum likelihood methods of estimation. The means and MSEs for the different estimators will be reported from these experiments. We can see from Table 1 when n is increase MSE is decrease.

Applications

This section provides an application to show how the MBGz distribution can be applied in practice. We compare MBGz to Exponentaited Generalized Weibull-Gompertz distribution (EGWGz) by [START_REF] El-Bassiouny | Exponentiated Generalized Weibull-Gompertz Distribution with Application in Survival Analysis[END_REF] and other well known distributions in literature, Kumaraswamy-Gompertz (Kw-Gz), beta Gompertz(BGz) and Gompertz (Gz) models. The MLEs are computed using Quasi-Newton Code for Bound Constrained Optimization and the log-likelihood function evaluated. The goodness-of-fit measures, Anderson-Darling (A*), Cramer-von Mises (W*), Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and log-likelihood ( ˆ ) values are computed. The lower the values of these criteria, the better the fit. The value for the Kolmogorov Smirnov (KS) statistic and its p-value are also provided. The required computations are carried out in the R software.

Data set 1

The first data set represents the time to failure(103 h) of turbocharger of one type of engine given in [START_REF] Xu | Application of Neural Networks in forecasting Engine Systems Reliability[END_REF]. The data set is as follows: 0.0312, 0.314, 0.479, 0.552, 0.700, 0.803, 0.861, 0.865, 0.944, 0.958, 0.966, 0.977, 1.006, 1.021, 1.027, 1.055, 1.063, 1.098, 1.140, 1.179, 1.224, 1.240, 1.253, 1.270, 1.272, 1.274, 1.301, 1.301, 1.359, 1.382, 1.382, 1.426, 1.434, 1.435, 1.478, 1.490, 1.511, 1.514, 1.535, 1.554, 1.566, 1.570, 1.586, 1.629, 1.633, 1.642, 1.648, 1.684, 1.697, 1.726, 1.770, 1.773, 1.800, 1.809, 1.818, 1.821, 1.848, 1.880, 1.954, 2.012, 2.067, 2.084, 2.090, 2.096, 2.128, 2.233, 2.433, 2.585, 2.585

Data set 2

The second data set was given by [START_REF] Badar | Statistical aspects of fiber and bundle strength in hybrid composites[END_REF]. It corresponds to a single fiber with 20 and 101 mm of gauge length, respectively. The data set is as follows: 1.6, 2.0, 2.6, 3.0, 3.5, 3.9, 4.5, 4.6, 4.8, 5.0, 5.1, 5.3, 5.4, 5.6, 5.8, 6.0, 6.0, 6.1, 6.3, 6.5, 6.5, 6.7, 7.0, 7.1, 7.3, 7.3, 7.3, 7.7, 7.7, 7.8, 7.9, 8.0, 8.1, 8.3, 8.4, 8.4, 8.5, 8.7, 8.8, 9.0 Tables 2 and4 list the maximum likelihood estimates (and the corresponding standard errors in parentheses) of the unknown parameters of the MBGz distribution for Data Set 1 and Data Set 2 respectively. Tables 3 and5 show the statistics AIC, BIC, W*, A*, KS, P-Value values for all the considered models. We then see that the proposed MBGz model fits these data better than the other models. The MBGz model may be an interesting alternative to other models available in the literature for modeling positive real data. To complete this fact, PP, QQ, epdf and ecdf plots of the MBGz distribution given in Figures 2 and3 for Data Set 1 and Data Set 2 respectively.

  8 α = 0.1 a = 2 b = 2 c = 0.8 λ = 0.9 α = 1 a = 0.1 b = 0.5 c = 0.01 λ = 1 α = 0.8 a = 3 b = 0.5 c = 0.5 λ = 1 α = 1 a = 0.2 b = 07 α = 0.1 a = 3 b = 5 c = 0.1

Figure 2 :Figure 3 :

 23 Figure 2: PP, QQ, epdf and ecdf plots of the MBGz distribution for Data Set 1.

Table 1 :

 1 The MLEs and MSEs of MBGz distribution.

	n	Parameters Initial	MLE	MSE	Initial	MLE	MSE
	50	a	3.0	3.0024 0.5057	2.5	2.6424 0.1737
		b	1.5	1.6409 0.1499	1.5	1.5219 0.0400
		c	0.5	0.4941 0.0008	0.5	0.5050 0.0004
		α	0.5	0.5422 0.0198	0.5	0.5291 0.0116
		λ	0.5	0.5241 0.0387	0.5	0.5235 0.0122
	100	a	3.0	3.0778 0.2458	2.5	2.5060 0.0754
		b	1.5	1.6083 0.0779	1.5	1.5373 0.0291
		c	0.5	0.4986 0.0003	0.5	0.4991 0.0003
		α	0.5	0.5572 0.0147	0.5	0.5123 0.0029
		λ	0.5	0.4926 0.0126	0.5	0.5035 0.0070
	150	a	3.0	2.9041 0.1015	2.5	2.5125 0.0284
		b	1.5	1.6159 0.0485	1.5	1.5232 0.0088
		c	0.5	0.4940 0.0002	0.5	0.5002 0.0001
		α	0.5	0.5477 0.0094	0.5	0.5137 0.0015
		λ	0.5	0.4694 0.0072	0.5	0.4968 0.0015
	50	a	1.5	1.4706 0.0325	1.5	1.5435 0.0641
		b	1.8	1.7764 0.0639	1.8	1.7838 0.0955
		c	0.5	0.5054 0.0013	1.5	1.5285 0.0203
		α	0.5	0.4833 0.0029	0.5	0.4895 0.0008
		λ	0.5	0.5488 0.0160	0.5	0.5364 0.0118
	100	a	1.5	1.5138 0.0201	1.5	1.5194 0.0224
		b	1.8	1.8177 0.0380	1.8	1.8309 0.0451
		c	0.5	0.5004 0.0007	1.5	1.5010 0.0047
		α	0.5	0.5007 0.0023	0.5	0.5011 0.0005
		λ	0.5	0.5106 0.0059	0.5	0.5036 0.0028
	150	a	1.5	1.5313 0.0102	1.5	1.4690 0.0094
		b	1.8	1.8152 0.0194	1.8	1.8396 0.0258
		c	0.5	0.5055 0.0003	1.5	1.4864 0.0017
		α	0.5	0.5173 0.0022	0.5	0.5007 0.0004
		λ	0.5	0.5044 0.0034	0.5	0.4943 0.0009

Table 2 :

 2 MLEs (standard errors in parentheses) for Data Set 1.

	Distribution			Estimates
	MBGz(λ, α, a, b, c)	0.0085	2.5537	1.0737	1.3153	5.0687
		(0.0067) (0.5727)	(0.3197)	(0.8933) (3.3003)
	EGWGz(λ, a, b, c, β)	3.2078	2.4598	0.0203	1.8974	0.5460
		(1.2099) (0.6498)	(0.0531)	(1.8193) (0.2430)
	KwGz(a, b, c, d, θ)	0.1861	1.4948	1.4909	0.9811
		(0.3130) (0.5076)	(0.4735)	(2.4368)
	BGz(a, b, θ, v)	0.3144	1.5591	1.4798	0.4966
		(0.4283) (0.3658)	(0.4543)	(0.8692)
	Gz(λ, α)	0.0841	1.8811
		(0.0268) (0.2043)

Table 3 :

 3 The AIC, BIC, W*, A*, KS, P-Value values for Data Set 1. Value MBGz 50.0387 110.0776 118.2481 0.0328 0.2745 0.0539 0.9889 EGWGz 52.6888 115.3776 126.5482 0.0706 0.5341 0.0785 0.7885 KwGz 51.2042 110.4084 119.3448 0.0529 0.4125 0.0640 0.9396 BGz 51.1518 110.3026 119.2399 0.0518 0.4057 0.0627 0.9484 Gz 53.9686 111.9374 122.4056 0.0819 0.5921 0.0810 0.7547

	Dist	ˆ	AIC	BIC	W*	A*	KS	P-

Table 4 :

 4 MLEs (standard errors in parentheses) for Data Set 2.

	Distribution			Estimates
	MBGz(λ, α, a, b, c)	0.0098	0.5270	0.8768	4.5635	0.1561
		(0.0116) (0.1599)	(0.3893)	(0.8862) (0.2442)
	EGWGz(λ, a, b, c, β)	0.0101	0.6077	0.1078	1.6929	0.6613
		(0.0141) (0.1506)	(0.3427)	(1.2539) (0.3379)
	KwGz(a, b, θ, v)	0.0133	0.2923	2.0164	13.7085
		(0.0120) (0.1641)	(0.7880)	(7.0208)
	BGz(a, b, c, d, θ)	0.0125	0.1856	3.7622	2.0116
		(0.0100) (0.1601)	(2.5635)	(3.3802)
	Gz(λ, α)	0.0074)	(0.6243)
		(0.0035) (0.0748)

Table 5 :

 5 The AIC, BIC, W*, A*, KS, P-Value values for Data Set 2.

	Dist	ˆ	AIC	BIC	W*	A*	KS	P-Value
	MBGz 78.2184 168.1770 176.0214 0.0222 0.1840 0.0707 0.9888
	EGWGz 79.3744 168.5489 178.5933 0.0479 0.2922 0.0821 0.9623
	KwGz 80.7197 169.4395 176.1950 0.0430 0.3326 0.0966 0.8489
	BGz	82.9924 173.9849 180.7404 0.0922 0.6736 0.1080 0.7389
	Gz	80.9566 168.9234 177.2911 0.0359 0.2335 0.0903 0.8299

Empirical and theoretical dens.


A I Buzdin

L N Bulaevskit

Ferromagnetic film on the surface of a superconductor: possible onset of inhomogeneous magnetic ordering

A system consisting of a bulky superconductor and a thin ferromagnetic metallic film on its surface is considered. It is shown that under certain conditions the magnetic ordering in the film takes the form of a small-scale domain structure. This ordering corresponds to a minimum system energy if account is taken of the interaction between the magnetic moments and the superconducting electrons that penetrate into the film from the bulky superconductor.

I.INTRODUCTION

The strong exchange field (exceeding as a rule by sever al orders the superconducting critical temperature Tc) pres ent in a ferromagnet suppresses superconductivity via the paramagnetic effect. 1 For the same reason, the supercon ducting order parameter is radically decreased near the con tact of a superconductor with a ferromagnetic metal. The principal role is played here not by the proximity itself to the normal metaV but precisely by the presence of a strong ex change field in the ferromagnet and by the electron penetra tion, due to the proximity effect, into this region of the strong exchange field.

The influence of the exchange field on the superconduc tivity is greatly weakened when a small-scale (compared with the superconducting correlation length So= 0.18 Vp/Tc) domain structure is produced in the ferromagnet the exchange field is effectively averaged over Cooper-pair dimensions. Precisely such domain coexistence phases were predicted for ferromagnetic superconductors, and the ex perimental data agreed with this prediction. 3 For this reason, inhomogeneous magnetic ordering, rather than ferromagnetism, should likewise occur in a thin ferromagnetic-metal film of thickness d (on the order of sev eral dozen Angstroms) sputtered on the surface of a super conductor, provided that the Curie temperature 0 < Tc. In the case Q� s 0-1, where Q is the characteristic wave vector of such a magnetic structure, the oscillating magnetic field in the film suppresses less the superconductivity at the inter face. The gain in the superconducting-condensation energy near the interface compensates here for the magnetic-energy loss due to formation of an inhomogeneous magnetic struc ture in the film.

A more interesting situation can be expected if 0 > Tc.

In this case, below the temperature Tc, in view of the weak ness of the superconducting pairing, the film should be in a It is interesting to note that a suitable object for the study of the transition considered may be vanadium. Onset of surface ferromagnetism in vanadium was observed recent ly.4 At low temperatures vanadium becomes superconduct ing and one can expect a domain magnetic structure rather than ferromagnetism to be produced on the surface.

The considered transition in a ferromagnetic film on the surface of a superconductor can be observed by magneto optic methods or with the aid of tunnel measurements. The use of a ferromagnetic film as the weak-coupling link of a Josephson structure would make feasible a Josephson junc tion with variable coupling-the transparency should in crease jumpwise when the film goes over into the domain phase.

BASIC EQUATIONS

We consider the case when a superconductor surface is coated with a thin ferromagnetic-metal film of thickness d that is small compared with the superconducting correlation length s ( T = 0), and assume that the easy axis of the ferro magnet is parallel to the plane of the film (if the easy axis is perpendicular to the film plane a small-scale domain struc ture exists even in the absence of superconductivity). We assume that the film can have two magnetic states-ferro magnetic without domains (or with domains very large compared with so>. and one with small-scale domain struc ture, see Fig. 1. We designate these states by F and DS.

The details of the F-DS transition in a ferromagnetic film on the surface of a superconductor depend on many parameters of the system: on the homogeneity of the film, on the conditions of the contact on the ferromagnet-supercon ductor interface, on the density of the electronic states in the film and in the bulk, and others. To illustrate the physics of the phenomenon, we confine ourselves therefore to a patent ly greatly simplified model. We assume that the dimension less Cooper-pairing constant A. is the same in the film and in the bulk, and the electron mean free path in the film is I < So •

We confine ourselves also to a situation in which the transi tion to the DS phase takes place near the superconducting transition. In the context of the BCS theory, the initial Ham iltonian of the system considered can be written in the form The solution of the equation that follows from ( 3) for a(x) is well known:

(4)

where $2 = 1]14m!r!, and the constant c is determined by the condition (2), viz, sinh(2C) = 2 1 1 2 /y$ (we have taken it into account that d <to$).

The influence of the magnetic structure of the film on the superconductivity is in fact determined in this case only by the parameter y. For the ferromagnetic state Fin the film, the superconducting order parameter (under the condition hdsolvpl> 1) is zero on the boundary, a( d) = 0, i.e., r-00.

The solution ( 4) for the F phase is

L1 (x) =L10 th(x/2'1'£).
The condensation-energy loss (per unit surface area) due to the suppression of the superconductivity amounts in this case to5

"t '

2'/, 6F=n-� � 3 .

(5)

If the DS phase is produced, the superconductivity of the interface is not so strongly suppressed (in view of the averaging of the exchange field), and the properties of the DS phase (period, exchange-field amplitude) are deter mined just by the parameter y in ( 2). The condensation energy loss is decreased and, as can be easily seen, equals A contribution to the total energy of a system with aDS phase is made also by the inhomogeneity of the magnetic structure in the film: the presence of a domain wall entails an energy loss Ew = 0nd(a/b), where b is the domain thick ness (half the DS-structure period), a is the magnetic corre lation length and is of the order of the interatomic distance, and 0 is of the order of the Curie temperature (the electron and magnetic-atom densities are assumed to be approxi mately equal). The expression forE w indicates that the en ergy loss per magnetic atom in a domain wall of approxi mately atomic thickness is � 0, i.e., the magnetic anisotropy is regarded as strong-in the case of weak magnetic anisot ropy the domain-wall energy is somewhat lower and a gener alization to include this case is easy.

The total system energy in the presence of aDS phase in the film can be written in the form

BOUNDARY CONDITION

Assuming the mean free path of the electrons in the ferromagnetic film to be small, we use the Usadel equations8

to describe the superconductivity in the DS-phase region ( 0 < x < d, see Fig. 1 ) .

In an approximation linear in the function F (the anom-

  ferromagnetic phase and suppress strongly the supercon ductivity in a region of orders( T) near the boundary. With further decrease of the temperature, the superconducting condensation energy increases, as does also the energy loss due to the suppression of the superconductivity near the boundary. It may be expedient as a result to break up the films into minute domains: energy loss due to domain-wall formation is offset by the superconducting-condensation en ergy gain near the boundary (the domain phase suppresses the superconductivity more weakly). The temperature T vs at which the film goes over into the domain phase decreases with increase of the film thickness and with decrease of the superconducting correlation length So of the superconduc tor. The most suitable for observation of the transition are pure superconductors with large So• and also superconduc tors with high condensation energy (high Tco high elec tronic-state density).

  FIG. I. de= J d3r{1jl+ (r)L ljl(r) + ,1 (r)ljl+ (r) iay \jl+ (r) 2m -,1'(r)1jJ(r)iay1jJ(r)+1jJ'(r)h(r)mjJ(r)+/L1(r) 1'/A. }+de,c, (1) where a(r) is the superconducting order parameter, 1/J(r) is a spin or, er are Pauli matrices, h ( r) = h is the exchange field acting on the electrons in the film ( h ( r) = 0 inside a bulky superconductor), and the term :lt' se takes into account elec tron scattering by nonmagnetic impurities. The influence of the magnetic induction on the superconductivity is shown in Ref. 3 to be negligibly small under ordinary conditions com pared with the exchange-field influence, and will be neglect ed here. To describe the behavior of the order parameter in a superconductor, except in region near the film, we can use in lieu of ( 1) the Ginzburg-Landau function Y (see, e.g., Ref.5), supplemented by the boundary condition a' (x) I a (x)-x�d = 'Y•

  the energy of the ferromagnetic ordering in the film, and E5 is the superconducting-condensation energy in the absence of a magnetic film on the surface. The values of Em and E5 are the same for both the F and DS phases, are independent b, and their specific form is immaterial. The total energy for the ferromagnetic state in the film is The period of the domain structure is determined from the condition a E DS I ab = 0, and the condition for the transition from the ferromagnetic to the DS structure is the equality Evs = Ep. The parameter C in (7) is determined by the boundary condition (2), where r depends on b. A complete solution of the problem calls thus for find ing the boundary condition, i.e., for calculating the param eter r as a function of b.

alous Green's function integrated over the energies and ve locity directions on the Fermi surface) , the Usadel equation in the presence of the exchange field h is '[ro+ih(x, y)]F-D-1DV2F=tl(x, y), (9) where {J) = TTT ( 2n + 1 ) , D = vFI 13 is the diffusion coeffi cient , and the approximation linear in F is used because we are considering the temperature region near T c .

Rapid oscillations of the exchange field in the DS phase along they coordinate cause a weak dependence of ll. on y , which can be neglected , i.e. , tl:::: : ll.(x). Using for F( x , y) the representation F(x, y) = F,(x) + L, Fk(x)e;k", ( 10) where k=kn = 1rnlb, n = 1,2, ... , and using the condition vFilk 2 )-h)>. T c

, we obtain for F 0 (x) the equation

where

( 12)

As seen from ( 11) the action of aDS structure on the super conductivity is analogous to that of magnetic scattering with a reciprocal time Ts 1 given by Eq. ( 12). This result is gen eral for dirty superconductors , see Ref. 3. From Eq. ( 1) and from the self-consistency condition we find that the behavior of the superconducting order pa rameter in the film is described by the equation

here 1/J(x) is a digamma function and it follows from ( 3) that the solution of interest to us , which satisfies the condition /l. ' (x) = 0 on the interface with the vacuum , is of the form tl(x)=chqx, q2=21't,D.

The boundary condition ( 2) takes thus the form tl' (x)

We determine the equilibrim energy in the case of aDS phase in the film. From the condition aEDslab = 0 , taking into account in ( 7) the dependence of C on the period b of the structure [ Eqs. ( 12) and ( 15) ] , we get

Using this condition , we find that at the transition point

and since sinh 2C = 2 112 lrs. we have

Substituting this value of bin ( 16) , we get the relative transi tion temperature v l 2'1• ( 2; ) '" th2C 't o_ ,, ,= d

It is taken into account here that s = s0r 0 1 1 2 • In order of magnitude , as follows from ( 18) , we have

Our analysis is qualitatively valid also at low temperatures , therefore the general condition for the transition to the DS phase

For a ferromagnet with RKKY interaction we have ®:::: :: h ) , and the destructive action of the exchange field in the DS phase turns out to be relatively weak ( Ts

). Note that we have disregarded magnetic scat tering in the film. Below the Curie point in the presence of easy-axis magnetic anisotropy , magnetic scattering is strongly suppressed. In principle , however , its influence can be easily taken into account , and furthermore the foregoing conclusions remain qualitatively in force even for a recipro cal magnetic scattering timer;;; 1-®.

It follows from our analysis that the most favorable conditions for the observation of the DS -F transition are produced in superconductors with a large correlation length So• Furthermore , the effect is enhanced by the high density of the electronic states in the superconductor.

We have considered a system in which the magnetic film is on the surface of a superconductor. A similar situa tion can arise also in a superlattice consisting of thick super conductor slabs separated by thin ferromagnetic-metal lay ers. The F-DS transition in the layer leads here to a first-order transition from two-dimensional to three-dimen sional superconductor behavior.


Nicolas Aragon
email: nicolas.aragon@unilim.fr

Philippe Gaborit
email: philippe.gaborit@unilim.fr

Adrien Hauteville
email: adrien.hauteville@unilim.fr

Jean-Pierre Tillich
email: jean-pierre.tillich@inria.fr

A New Algorithm for Solving the Rank Syndrome Decoding Problem

Abstract-In this paper, we propose an improvement of the attack on the Rank Syndrome Decoding (RSD) problem found in [START_REF] Gaborit | On the complexity of the rank syndrome decoding problem[END_REF], usually the best attack considered for evaluating the security of rank based cryptosystems. For H a full-rank (n -k) × n matrix over Fqm and e ∈ F n q m of small norm r, the RSD problem consists in recovering e from s = He T . In our case, the norm of a vector over Fqm is defined by the dimension of the Fq-subspace generated by its coordinates. This problem is very similar to the Syndrome Decoding problem in the Hamming metric (only the metric and the field of the coefficients are different) and the security of several cryptosystems relies on its hardness, like McEliece-based PKE [START_REF] Gaborit | New results for rankbased cryptography[END_REF], [START_REF] Loidreau | A new rank metric codes based encryption scheme[END_REF] or IBE [START_REF] Gaborit | Identity-based encryption from rank metric[END_REF]. Our attack is in O (nk) 3 ,k+1

[1], [START_REF] Ourivski | New technique for decoding codes in the rank metric and its cryptography applications[END_REF]. In particular in the case m ≤ n, our attack permits to obtain an exponential gain in q m(1-R) for R = k/n the rate of the code.

We give examples of broken parameters for recently proposed cryptosystems based on LRPC codes or Gabidulin codes. Our attack does not fully break these cryptosystems but implies larger parameters for the same security levels.

I. INTRODUCTION

The hardness of the problem of decoding a linear code makes its use very attractive for post-quantum cryptography. Indeed it has been proven to be NP-complete for the Hamming metric [START_REF] Berlekamp | On the inherent intractability of certain coding problems[END_REF]. One of the first cryptosystem partly based on this problem is the McEliece cryptosystem, whose instances using the family of Goppa codes or the family of MDPC codes still remain unbroken today. One of the main drawback of these cryptosystems are the relatively large public key sizes. A solution to decrease their size is to change the metric used for the code.

The rank metric F q m -linear codes has been introduced in [START_REF] Gabidulin | Theory of codes with maximum rank distance[END_REF] and their first use for the cryptosystem of McEliece was in 1991 [START_REF] Gabidulin | Ideals over a non-commutative ring and their applications to cryptography[END_REF]. This cryptosystem used the family of Gabidulin codes and a structural attack has been found by Overbeck in [START_REF] Overbeck | A new structural attack for GPT and variants[END_REF]. However another family, the Low Rank Parity Check (LRPC) codes, was proposed for the McEliece cryptosystem in [START_REF] Gaborit | New results for rankbased cryptography[END_REF]. These codes are similar to the MDPC codes in the Hamming metric [START_REF] Misoczki | MDPC-McEliece: New McEliece variants from moderate density parity-check codes[END_REF] or the NTRU cryptosystem in the euclidean metric [START_REF] Hoffstein | NTRU: A ring-based public key cryptosystem[END_REF] and have a poor algebraic structure, which ensures a good resilience against structural attacks.

Since the security of the McEliece cryptosystem relies directly on the hardness of decoding a random linear code, it is very important to evaluate the complexity of generic algorithms to solve it. One of the first algorithm [START_REF] Ourivski | New technique for decoding codes in the rank metric and its cryptography applications[END_REF] has a complexity of O (wm) 3 q (w-1)(k+1)+2 for an error of weight w in a code of length n and of dimension k over the field F q m . As we can see, this complexity does not take into account the length of the code. The GRS algorithm [START_REF] Gaborit | On the complexity of the rank syndrome decoding problem[END_REF] gives a complexity of O (n -k) 3 m 3 q (w-1) (k+1)m n in the case n m and this article introduces another algorithm based on the qpolynomials of complexity O (wm) 3 q r w (k+1)(w+1)-n-1 w . In this paper, we improve the complexity of the GRS algorithm and we obtain a complexity of O (n-k) 3 m 3 q w (k+1)m n -m .

II. BACKGROUND ON THE RANK METRIC

A. Definitions

Let us begin with the definitions of matrix codes and rank metric.

Definition II.1 (Linear matrix codes). A linear matrix code C of length m × n and dimension K over F q is a subspace of dimension K of M m×n (F q ). We say C is an [m × n, K] q linear matrix code, or simply an [m × n, K] code if there is no ambiguity.

We define the rank metric over matrix codes as such:

• the distance between two words A and

B is d R (A, B) def = Rank(A -B). • the weight of a word A is |A| r def = Rank(A) = d R (A, 0).
An important family of matrix code are the F q m -linear codes.

Definition II.2 (F q m -linear codes). A F q m -linear code C of length n and dimension k is a subspace of dimension k of F n q m . We say C is an [n, k] q m linear code, or simply an [n, k] code if there is no ambiguity.

A natural way to define the rank metric over F q m -linear codes is to consider the matrix associated to a word of F n q m . Definition II.3 (Associated matrix). Let x = (x 1 , . . . , x n ) ∈ F n q m and β = (β 1 , dots, β m ) a basis of F q m over F q .

x can be represented by an m × n matrix M x such that its i th column represents the coordinates of x i in the basis β.

x

↔    x 11 . . . x 1n . . . . . . x m1 . . . x mn    with x j = m i=1 x ij β i for all j ∈ [1..n]. By definition, |x| r def = |M x | r and d R (x, y) = d R (M x , M y ).
These definitions do not depend on the choice of the basis.

Definition II.4 (Matrix code associated to an F q m -linear code). Let C be an [n, k] q m linear code. The matrix code C M associated to C (with respect to a basis of F q m over F q ) is the code

{M x : x ∈ C}.
It is an [m × n, km] q matrix code.

An [n, k] q m linear code C can be represented by a generator matrix G ∈ F k×n q m where each row g i of G is an element of a basis of C or by a parity-check matrix

H ∈ F (n-k)×n q m such that c ∈ C ⇐⇒ cH T = 0.
In the first case, C is viewed as the image of a morphism and in the second case as the kernel of a morphism. If G (respectively H) is of the form (I k |A) (respectively (B|I n-k )), we say that it is under systematic form.

The advantage of F q m -linear codes with respect to matrix codes of the same parameters is that they have a much more compact representation. Indeed, under their systematic form, an [n, k] q m F q m -linear can be represented by

(n -k)k coefficients in F q m which is (n -k)km coefficients in F q whereas an [m × n, km] q matrix code can be represented by (nm -km)km = (n -k)m 2 coefficients in F q .
Definition II.5 (Support of a word). Let x ∈ F n q m . The support of x is the F q -subspace of F q m generated by the coordinates of x. It is denoted Supp(x).

Supp(x) = x 1 , . . . , x n Fq
The weight of a word is equal to the dimension of its support.

The number of supports of dimension w in F q m is denoted by the Gaussian coefficient m w q and is equal to:

m w q = w-1 i=0 q m -q i q w -q i
In practice Gaussian coefficients are approximated by q w(m-w) . In the case of matrix codes, we can consider the subspace of F m q generated by the columns of the matrix, called the columns support, or the subspace of F n q generated by its rows, called the row support.

B. Hard problem in rank metric

Rank-based cryptography relies on difficult problems in rank metric. These problems are the same as in the Hamming metric, but with the rank metric replacing the Hamming metric. In this subsection, we only consider F q m -linear codes but all the problems are defined exactly the same way for linear matrix codes.

The first one corresponds to the decoding problem by syndromes.

Definition II.6 (Rank Syndrome Decoding (RSD) problem).

Let H be the parity-check matrix of an

[n, k] code, s ∈ F n-k q m
and w an integer. The RSD problem consists in finding e ∈ F n q m such that:

He T = s T |e| r = w (1)
The second problem is the search for small-weight codewords.

Definition II.7 (Small-weight codeword problem). Let C be an [n, k] code and w an integer. The problem consists in finding a codeword c ∈ C such that |c| r = w.

Remark: If H is a parity-check matrix of C, the smallweight codeword problem corresponds to an instance of the RSD problem with s = 0.

III. DESCRIPTION OF THE ATTACKS ON THE RSD PROBLEM

Our algorithm is an improvement of the GRS algorithm [START_REF] Gaborit | On the complexity of the rank syndrome decoding problem[END_REF]. We first recall the description of this algorithm in the first subsection then we present our improvement in the second one.

A. The GRS algorithm

The general idea to solve the RSD problem is to find a subspace F such that Supp(e) ⊂ F . Then we can express the coordinates of x in a basis of F and solve the linear system given by the parity-check equations. There are two possible cases we will describe more precisely.

• first case: n m.

Let F be a subspace of F q m of dimension r and (F 1 , . . . , F r ) a basis of F . Let us assume that Supp(e) ⊂ F .

⇒ ∀i ∈ [1..n], e i = r j=1 λ ij F j
This gives us nr unknowns over F q .

We can now rewrite the parity-check equations in these unknowns:

He T = s (2) ⇔      H 11 e 1 + • • • + H 1n e n = s 1 . . . . . . . . . H n-k,1 e 1 + • • • + H n-k,n e n = s n-k ⇔      n l=1 H 1l r j=1 λ lj F j = s 1 . . . . . . n l=1 H n-k,l r j=1 λ lj F j = s n-k (3)
Now, we need to embed these n -k equations over F q m into (n -k)m equations over F q . Let ϕ i the i th canonical projection from F q m on F q :

ϕ i : F q m → F q m i=1 x i β i → x i
We have:

He T = s ⇔ ∀i ∈ [1..m],      n l=1 r j=1 λ 1j ϕ i (H 1l F j ) = ϕ i (s 1 ) . . . . . . n l=1 r j=1 λ 1j ϕ i (H n-k,l F j ) = ϕ i (s n-k ) (4)
Since we assume Supp(e) ⊂ F , this system has at least one solution. We want more equations than unknowns in order to have only one solution with overwhelming probability. So we have the condition:

(n -k)m nr ⇔ r m -km n

If the system has a solution, we check if its rank is equal to w. Otherwise, we choose another F and restart at the beginning.

The average complexity of the algorithm is

O (n-k) 3 m 3 p
where p is the probability that Supp(e) ⊂ F . The probability p is equal to the number of subspaces of dimension w in a subspace of dimension r divided by the number of subspaces of dimension w in

F q m . p = r w q m w q ≈ q -w(m-r)
By taking r = m -km n we obtain a complexity of O (n -k) 3 m 3 q w km n • second case: m > n. In this case we consider the matrix code associated to the F q m -linear code. Instead of searching for a subspace which contains the columns support of the matrix of the error, we search for a subspace F of F n q which contains the rows support of the error. The rest of the algorithm is the same, the only differences are:

the number of unknowns is mr, which implies r n -k.

probability to choose F correctly is

  r w   q   n w   q ≈ q -w(n-r) .
Thus, the complexity in this case is O (n -k) 3 m 3 q wk . This algorithm does not use the F q m -linearity of the code. To exploit this structure, the main idea is to consider the code C = C + F q m e, where C is the code with parity-check matrix H.

The problem is reduced to the search for a codeword of weight w in C . If e is the only solution of the system 1, then the only codewords of C of weight w are of the form αe, α ∈ F * q m . These codewords are solutions of the system

H e T = 0 |e | r = w ( 5
)
where H is a parity-check matrix of C .

Thus, instead of looking for the support E of e, we can look for any multiple αE of the support. In [START_REF] Gaborit | On the complexity of the rank syndrome decoding problem[END_REF], the authors specialize one element of the support by searching for small weight codeword c such that 1 ∈ Supp(c) = E (but one can choose any non zero element of F * q m ). One can use this information to improve the probability that F ⊃ Supp(c) by choosing F such that 1 ∈ F . Let ϕ be the projection from the F q -subspaces of F q m which contains 1 to the subspaces of the quotient-space F q m /F q such that ϕ(V ) = V /F q . We have dim ϕ

(F ) = dim F -1 and dim ϕ(E) = dim E -1. Moreover, E ⊂ F if and only if ϕ(E) ⊂ ϕ(F ).
Thus, the probability is equal to

  w -1 r -1   q   w -1 m -1   q ≈ q (w-1)(m-r) .
We use the same techniques as before to obtain the system

     n l=1 r j=1 λ 1j ϕ i (H 1l F j ) = 0 . . . . . . n l=1 r j=1 λ 1j ϕ i (H n-k-1,l F j ) = 0 (6)
We take r = m(n-k-1) n = m -m(k+1) n in order to have more equations than unknowns. If F ⊃ E, then the kernel of system 6 is of positive dimension and we can compute a non-zero vector e solution of H e T = 0. Then if it is of rank w, we solve the equation He T = αHe T = αs T of unknown α to recover e. Finally, we get an average complexity of O (n -k) 3 m 3 q (w-1) (k+1)m n .

In the case m > n, we cannot use the attack of Ourivski and Johansson to obtain an average complexity of O (nk) 3 m 3 q (w-1)(k+1) [START_REF] Ourivski | New technique for decoding codes in the rank metric and its cryptography applications[END_REF].

B. Improvement of this algorithm

The idea of our improvement is still to search for a codeword of weight w in the code C . However, instead of choosing F such that 1 ∈ F , we choose a random F . If it contains a multiple αE of E, then we can compute the codeword αe of C as before. The probability that F ⊃ αE depends on the number of different subspaces of the form αE, α ∈ F * q m . Proposition III.1. Let E be an F q -subspace of F q m of dimension w. There are at most q m -1 q-1 subspaces of the form αE, α ∈ F * q m , with equality if w and m are coprime. Proof. Let S be the cardinality of the set of the subspaces of this form. Let α and β ∈ F q m such that αβ -1 ∈ F * q By linearity, it is obvious that αE = βE, hence S q m -1 q-1 . To prove the equality, let α ∈ F * q m be such that E = αE. By induction, we have F q (α)E = E thus E is an F q (α)-subspace of F q m . Let m = [F q (α) : F q ]. We have w = dim Fq E = m dim Fq(α) E hence m divides w and m. If GCD(m, w) = 1 then α ∈ F * q , which concludes the proof. This also proves that if the number S of different subspaces of the form αE is less than q m -1 q-1 then E is an F q m -subspace of dimension w m of F q m . The probability of such an event is negligible, so we can assume S = q m -1 q-1 .

The probability p that F of dimension r = m -m(k+1) n contains a subspace of the form αE, α ∈ F * q m can be approximated by the product of S by the probability that F contains a fixed subspace of dimension w. This approximation is correct if q m -1 q -1 r w q m w q ⇔ q m+w(r-w) q w(m-w) .

Notice that the bound q m -1 q-1 on S is an upper bound, the previous discussion shows that in practice the practical value of S can be approximated by this bound, which is used for evaluating the complexity of the attack. Now a value of S lower than the bound may permit to improve a little bit the attack but is unlikely to happen. The likely value of S associated with the probability that F contains contains a fixed subspace of dimension w leads us to the following theorem:

Theorem III.2. Let H ∈ F (n-k)×n q m , s ∈ F n-k q m
and w an integer. Our algorithm solves the RSD problem II.6 with average complexity of

O (n -k) 3 m 3 q w (k+1)m n -m operations in F q .
An implementation of the attack showed experimental values very close to the theoretic ones, which confirms our approximation. In the case m n, our algorithm is always better than [START_REF] Gaborit | On the complexity of the rank syndrome decoding problem[END_REF]. The gain in the exponent is equal to

m -(k+1)m n = m -mR ≈ m(1 -R )
where R is the rate of C . In the case m > n our algorithm may also be faster for some parameters.

IV. EXAMPLES OF BROKEN PARAMETERS

In this section, we present two public key encryption cryptosystems: the Loidreau cryptosystem [START_REF] Loidreau | A new rank metric codes based encryption scheme[END_REF] (based on Gabidulin codes masked by an LRPC matrix) and the LRPC cryptosystem [START_REF] Gaborit | New results for rankbased cryptography[END_REF]. Some of their parameters are broken by our attack. They are based on the McEliece cryptosystem. To encrypt a message m, one computes the cipher text c = mG pub + e where G pub is the public generator matrix of the masked code ans e is an error of wight w. To decrypt, the receiver removes the mask and decodes the error.

Under the assumption that the masked code is indistinguishable from a random code, an attacker has to solve an instance of the RSD problem of weight w to recover the message.

The following tables give some parameters which are broken by our attack:

• Parameters from Loidreau's cryptosystem [ For all these cryptosystems our attack improves on the best attack on which their security is based, the gain is a little less than m(1 -R) since the authors probably took a security margin on their parameters, but our attack still permits to break these parameters. Notice that some parameters of these cryptosystems are not affected by our attack since it is not necessarily the best attack if m > n.

Our attack also breaks the parameters of the IBE in [START_REF] Gaborit | Identity-based encryption from rank metric[END_REF] but they have been updated in the eprint version of the paper.

In the process of the standardization of postquantum cryptosystems engaged by the NIST, the cryptosystems LAKE, LOCKER, Ouroboros-R, RankSign and RQC are based on rank metric codes and take account on this attack for the choice of their security parameters. See https://csrc.nist.gov/Projects/Post-Quantum-Cryptography/Round-1-Submissions for the documentation on these submissions.

  in Fq whereas the previous best attacks are in O (n -k) 3 m 3 q (w-1) min (k+1)m n