Next Article in Journal
Multicollinearity and Linear Predictor Link Function Problems in Regression Modelling of Longitudinal Data
Previous Article in Journal
Quantile Dependence between Crude Oil Returns and Implied Volatility: Evidence from Parametric and Nonparametric Tests
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Scale Mixture of Maxwell-Boltzmann Distribution

by
Jaime S. Castillo
1,
Katherine P. Gaete
1,
Héctor A. Muñoz
1,
Diego I. Gallardo
2,*,
Marcelo Bourguignon
3,
Osvaldo Venegas
4 and
Héctor W. Gómez
1
1
Departamento de Matemáticas, Facultad de Ciencias Básicas, Universidad de Antofagasta, Antofagasta 1240000, Chile
2
Departamento de Matemática, Facultad de Ingeniería, Universidad de Atacama, Copiapó 7820436, Chile
3
Departamento de Estatística, Universidade Federal do Rio Grande do Norte, Natal 59078-970, RN, Brazil
4
Departamento de Ciencias Matemáticas y Físicas, Facultad de Ingeniería, Universidad Católica de Temuco, Temuco 4780000, Chile
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(3), 529; https://doi.org/10.3390/math11030529
Submission received: 18 November 2022 / Revised: 28 December 2022 / Accepted: 4 January 2023 / Published: 18 January 2023
(This article belongs to the Section Probability and Statistics)

Abstract

:
This paper presents a new distribution, the product of the mixture between Maxwell-Boltzmann and a particular case of the generalized gamma distributions. The resulting distribution, called the Scale Mixture Maxwell-Boltzmann, presents greater kurtosis than the recently introduced slash Maxwell-Boltzmann distribution. We obtained closed-form expressions for its probability density and cumulative distribution functions. We studied some of its properties and moments, as well as its skewness and kurtosis coefficients. Parameters were estimated by the moments and maximum likelihood methods, via the Expectation-Maximization algorithm for the latter case. A simulation study was performed to illustrate the parameter recovery. The results of an application to a real data set indicate that the new model performs very well in the presence of outliers compared with other alternatives in the literature.

1. Introduction

The Maxwell-Boltzmann (MB) distribution was introduced by Maxwell [1] to describe the distribution of speeds of molecules at thermal equilibrium and nowadays is widely applied in many fields such as statistical physics, statistical mechanics and accounting theory, among others. The MB distribution has been discussed in many works in the literature, for example, Tyagi and Bhattacharya [2] and Bekker and Roux [3].
Some recent extensions of the MB distribution are discussed, for example, in Sharma et al. [4], Vivekanand et al. [5], Iriarte et al. [6], Dey et al. [7], Sharma et al. [8] and Segovia et al. [9]. Product distributions or independent random variable quotients are of great interest; for example, Shakil et al. [10] studied the X Y and X / Y distribution, where X and Y are independent random variables that have MB and Rayleigh distributions respectively.
A random variable V follows the MB distribution with scale parameter β , denoted as V M B ( β ) , if its probability density function (pdf) and cumulative distribution function (cdf) are given by
f V ( v ; β ) = 4 β 3 / 2 π v 2 e β v 2 a n d F V ( v ; β ) = 2 π γ 3 2 , β v 2 ,
respectively, where v , β > 0 and γ ( a , v ) = 0 v t a 1 e t d t is the incomplete gamma function.
An extension of the MB distribution, called slash Maxwell-Boltzmann (SMB), was proposed by Acitas et al. [11]. The SMB distribution, denoted as S M B ( β , q ) , is defined as
Z = V U 1 / q , q > 0 ,
where V and U are independent random variables with MB( β ) and U ( 0 , 1 ) distributions, respectively, and q > 0 . As β in the MB distribution is a scale parameter, the motivation to define Z is the introduction of a shape parameter q, which makes this distribution more flexible. Note also that the model is parsimonious, being an alternative to traditional models such as the Weibull or gamma distributions, for example. The pdf of Z is
f z ( z ; β , q ) = 2 q Γ ( 1 / 2 ) β Γ q + 3 2 z β ( 1 + q ) G z 2 ; q + 3 2 ; β 2 , z , β , q > 0 ,
where Γ ( · ) denotes the gamma function and G ( · ; a ; b ) the cdf of a gamma distribution with shape and scale parameters a and b, respectively.
The main object of this article is to study an extension of the MB distribution with a greater range of the kurtosis coefficient, in order to use this new distribution to model datasets with atypical observations. We employ the slash methodology, as in the SMB version proposed by Acitas et al. [11]. Other authors have successfully applied the slash methodology. To name a few, Reyes et al. [12] obtained a generalization of the Birnbaum-Saunders (BS) distribution and Astorga et al. [13] introduced an extension on the power Muth distribution. In this work, we will show that the new distribution has heavier tails than the SMB distribution. Furthermore, this new distribution can be represented as a mixture of scales that allows us to perform simulation studies and obtain maximum likelihood (ML) estimators by means of the Expectation-Maximization (EM) algorithm.
The paper is organized as follows. Section 2 contains the representation of this model, and we generate the density of the new distribution. We present the scale mixture property, and the closed expressions for pdf, cdf, moments and coefficients of skewness and kurtosis, hazard and survival functions, and the Rényi entropy. Section 3 contains the inference, where we obtain the moments and ML estimators, and the implementation of the EM algorithm. In Section 4 we carry out a simulation study to assess the performance of the ML estimators in finite samples. Section 5 presents an application, comparing the fit of the scale mixture Maxwell–Boltzmann (SMMB) with the SMB, Weibull (W) and BS distributions to a real data set. Finally, Section 6 presents some conclusions.

2. Definition and Properties

Our proposal is based on the generalized gamma (GG) distribution introduced by Stacy [14]. The pdf for this model is given in Definition 1.
Definition 1.
A random variable Z follows the three-parameter GG distribution, denoted by Z G G ( a , d , p ) , if its pdf is
f Z ( z ; a , d , p ) = p a d Γ d p z d 1 e a z p ,
with a > 0 , d > 0 , p > 0 and z > 0 .
Definition 2.
A random variable X follows a SMMB distribution with scale parameter β > 0 and shape parameter q > 0 , denoted by X S M M B ( β , q ) , if X can be expressed as the ratio
X = V W ,
where V M B ( β ) and W G G ( 1 , q , 2 ) , both independent.
The following Proposition presents the pdf for the SMMB model.
Proposition 1.
Let X S M M B ( β , q ) with β > 0 and q > 0 . Then the pdf of X is
f X ( x ; β , q ) = 2 β 3 / 2 x 2 B q 2 , 3 2 1 + β x 2 q + 3 2 ,
where x > 0 and B ( · , · ) denotes the beta function.
Proof. 
Using the representation given in (2) and computing the Jacobian transformation, we have that:
X = V W Z = W V = X Z W = Z J = v x v z w x w z = z x 0 1 = z .
Then,
f X , Z ( x , z ) = | J | f V , W ( x z , z ) = 8 β 3 / 2 x 2 π Γ ( q / 2 ) z q + 2 exp z ( 1 + β x 2 ) , x > 0 , z > 0 .
The marginal pdf of X is:
f X ( x ; β , q ) = 8 β 3 / 2 x 2 π Γ ( q / 2 ) 0 z ( q + 3 ) 1 exp z ( 1 + β x 2 ) d z , x > 0 .
Note that the integrand in (3) is related to the pdf of a GG ( 1 + β x 2 ) 1 / 2 , q + 3 , 2 distribution. Therefore, the desired result follows.    □
The following proposition presents the cdf of the SMMB model, which is an expression involving the hypergeometric function that is defined by the power series:
2 F 1 ( a , b , c ; x ) = n = 0 ( a ) n ( b ) n ( c ) n x n n ! ,
where ( a ) n = 1 , n = 0 ; a ( a + 1 ) ( a + n 1 ) , n > 0 . For more details we refer the reader to Abramowitz and Stegun [15].
Proposition 2.
Let X S M M B ( β , q ) . Then the cdf of X is given by
F X ( x ; β , q ) = 2 β 3 / 2 x 3 3 B q 2 , 3 2 2 F 1 q + 3 2 , 3 2 , 5 2 ; β x 2 ,
where x > 0 , β > 0 and q > 0 .
Proof. 
We can write
F X ( x ; β , q ) = 0 x 2 β 3 / 2 t 2 B q 2 , 3 2 ( 1 + β t 2 ) q + 3 2 d t = 2 β 3 / 2 B q 2 , 3 2 0 x t 2 ( 1 + β t 2 ) q + 3 2 d t .
Hence, we make the following change of variable u = t 2 , to obtain
F X ( x ; β , q ) = β 3 / 2 B q 2 , 3 2 0 x 2 u 1 2 ( 1 + β u ) q + 3 2 d u .
Using the following result presented by Singh [16] (page 5)
0 t u a ( 1 + β u ) b d u = t a + 1 a + 1 2 F 1 b , a + 1 , a + 2 ; β t ,
in Equation (4), the result is shown.    □
Figure 1 shows the pdf and cdf for the SMMB ( β = 1 , q ) , for different values of q.
Remark 1.
For W G G ( 1 , q , 2 ) , it is possible to check that U = W 2 G ( q / 2 , 1 ) , i.e., the traditional gamma distribution with shape parameter q / 2 and rate 1. It therefore follows, from the properties of the inverse gamma model, that
E ( U 1 ) = 1 q / 2 1 a n d V ( U 1 ) = 1 ( q / 2 1 ) 2 ( q / 2 2 ) , i f q > 4 .
Provided that E ( U 1 ) 0 and E ( V 1 ) 0 , as q , it follows that
U 1 = 1 W 2 P 0 , a s q ,
where P denotes convergence in probability. Equivalently,
1 W P 0 , a s   q .
Based on the stochastic representation for the SMMB distribution in Equation (2), it follows that SMMB ( β , q ) P 0 , as q .

2.1. Lifetime Analysis

As the SMMB distribution is related to a non-negative and asymmetric variable, it can be used to model survival time data. In this section, the main features of interest in this field are studied. The survival and hazard functions for a SMMB model are provided in Corollaries 1 and 2.
Corollary 1.
Let X S M M B ( β , q ) . Then, the survival function of the X ( R X ) is given by
R X ( x ; β , q ) = 3 B q 2 , 3 2 2 β 3 / 2 x 3 2 F 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 3 B q 2 , 3 2 , x , β , q > 0 .
Corollary 2.
Let X S M M B ( β , q ) . Then, the hazard function of the X ( h X ) is given by
h X ( x ; β , q ) = 6 β 3 / 2 x 2 B q 2 , 3 2 ( 1 + β x 2 ) q + 3 2 3 B q 2 , 3 2 2 β 3 / 2 x 3 2 F 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 .
Figure 2 shows the survival and hazard functions. Additionally, we can see that the curve related to the hazard function is unimodal, and that as q grows, the curve has longer tails and extends over a greater range.
Note that considering η ( t ) : = f X ( t ) f X ( t ) = β ( q + 3 ) t 1 + β t 2 2 t for t > 0 , we have η ( t ) = β 2 ( q + 1 ) t 4 + β ( q + 7 ) t 2 + 2 t 2 ( 1 + β t 2 ) 2 . From this we obtain that
t 2 ( 1 + β t 2 ) 2 η ( t ) = β 2 ( q + 1 ) t 4 + β ( q + 7 ) t 2 + 2 ,
which implies that the zeros and signs of the η are same as those of the polynomial p ( t ) : = β 2 ( q + 1 ) t 4 + β ( q + 7 ) t 2 + 2 .
We note that t = q + 7 2 β ( q + 1 ) is a positive zero of the p ( t ) = 2 β t ( 2 β ( q + 1 ) t 2 + q + 7 ) , and p q + 7 2 β ( q + 1 ) = 12 β 2 ( q + 1 ) q + 7 2 β ( q + 1 ) + 2 β ( q + 7 ) = 4 β ( q + 7 ) < 0 . Therefore, at t = q + 7 2 β ( q + 1 ) a maximum is reached, with value: p q + 7 2 β ( q + 1 ) = ( q + 7 ) 2 4 ( q + 1 ) + 2 > 2 = p ( 0 ) . As t 0 = q + 7 + q 2 + 22 q + 57 2 β ( q + 1 ) is a zero of p ( t ) , we can conclude that
η ( t ) > 0 f o r t ( 0 , t 0 ) a n d η ( t ) < 0 f o r t ( t 0 , ) .
On the other hand, if we consider g ( t ) defined as in Glaser [17], in our case we have that
g ( t ) = ( β ( q + 1 ) t 2 2 ) ( 1 + β t 2 ) q + 1 2 t 3 π Γ ( q / 2 ) 4 β 3 / 2 Γ ( ( q + 3 ) / 2 ) 1 3 t 3 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 ) 1 ,
where 2 F 1 is the hypergeometric function. We note that g ( t ) tends to as t tends to 0 + ; and taking x = β t , g ( t ) it can be rewritten as
f q ( x ) = ( ( q + 1 ) x 2 2 ) ( 1 + x 2 ) q + 1 2 x 3 π Γ ( q / 2 ) 4 Γ ( ( q + 3 ) / 2 ) 1 3 x 3 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; x 2 ) 1 ,
so the number of positive zeros of g ( t ) is equal to that of f q ( x ) . In Table 1, the hypergeometric function 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 ) is shown for different values of q (see Weisstein [18]).
Based on g ( t ) graphs of the function defined in (5) for different q values shown in Figure 3, and using Theorem (d) part (i) in Glaser [17], we obtain that failure time density has an inverted bathtub shape.
The following proposition presents the scale mixture property for the SMMB distribution.
Proposition 3.
If X | V = v M B ( v 2 β ) and V G G ( 1 , q , 2 ) , is obtained X S M M B ( β , q ) .
Proof. 
We have that the joint function of ( X , V ) is f X , V ( x , v ) = f X | V ( x | v ) f V ( v ) , therefore, f X ( x ) is given by
f X ( x ; β , q ) = 0 f X , V ( x , v ) d v = 8 β 3 2 x 2 π Γ ( q 2 ) 0 v ( q + 3 ) 1 exp v ( 1 + β x 2 ) 1 2 2 d x .
Note that this function corresponds to the pdf of a random variable with GG ( ( 1 + β x 2 ) 1 2 , q + 3 , 2 ) distribution. Then the integral is equal to 1, consequently we have that
f X ( x ; β , q ) = 2 β 3 / 2 x 2 B q 2 , 3 2 1 + β x 2 q + 3 2 ,
where x > 0 and B is the beta function, also X S M M B ( β , q ) .    □
Remark 2.
Proposition 1 shows that the SMMB pdf has a closed expression, as occurs with its respective cdf given in Proposition 2. Proposition 3 shows that a SMMB distribution can also be obtained as a scale mixture of a MB and a GG distribution. This property is very important for generating random numbers and the implementation of the EM algorithm.

2.2. Moments

In this subsection, the following proposition shows the computation of the moments of a random variable with SMMB ( β , q ) distribution. Hence, it also displays the coefficients of skewness and kurtosis. For this, the following lemma will be useful.
Lemma 1.
Let W G G ( 1 , q , 2 ) with q > 0 . For r > 0 and r < q , then the r-th moment of W r is given by
E ( W r ) = Γ q r 2 Γ q 2 .
Proof. 
By definition,
E ( W r ) = 2 Γ ( q / 2 ) 0 w r q 1 e w 2 d w ,
and hence the result follows making the variable transformation u = w 2 .    □
In the next proposition, the moments of a SMMB model are given:
Proposition 4.
Let X S M M B ( β , q ) . Then the r-th moment of X is given by
μ r = E ( X r ) = 2 π β r / 2 Γ r + 3 2 Γ q r 2 Γ q 2 .
In particular,
μ 1 = 2 π β Γ q 1 2 Γ q 2 , q > 1 , μ 2 = 3 2 β Γ q 2 2 Γ q 2 , q > 2 ,
μ 3 = 4 π β 3 / 2 Γ q 3 2 Γ q 2 , q > 3 , μ 4 = 15 4 β 2 Γ q 4 2 Γ q 2 , q > 4 ,
V ( X ) = 3 π Γ q 2 2 Γ q 2 8 Γ q 2 2 2 2 π β Γ q 2 2 , q > 2 .
Proof. 
Using the stochastic representation for the distribution given in (2), we have that
E ( X r ) = E V W r = E ( V r ) E ( W r ) = 2 π β r / 2 Γ r + 3 2 Γ q r 2 Γ q 2 ,
where E ( V r ) = 2 π β r / 2 Γ r + 3 2 is the r-th moment of a MB ( β ) distribution and E ( W r ) was given in (6).    □
The asymmetry and kurtosis coefficients of the SMMB model are presented in the following proposition.
Proposition 5.
Let X S M M B ( β , q ) with β > 0 . Then the asymmetry and kurtosis coefficients of X are
β 1 = 2 2 4 π a 31 a 02 9 π a 12 a 21 a 01 + 16 a 13 27 π 3 a 23 a 03 216 π 2 a 22 a 02 a 12 + 576 π a 21 a 01 a 14 512 a 16 1 / 2 , f o r   q > 3 a n d β 2 = 15 π 2 a 41 a 03 128 π a 11 a 31 a 02 + 136 π a 12 a 21 a 01 192 a 14 9 π 2 a 22 a 02 48 π a 21 a 01 a 12 + 64 a 14 , f o r   q > 4 ,
respectively, where a i j = Γ q i 2 j .
Proof. 
Recall that by definition
β 1 = E [ ( X E ( X ) ) 3 ] ( V ( X ) ) 3 / 2 = μ 3 3 μ 1 μ 2 + 2 μ 1 3 μ 2 μ 1 2 3 / 2 , a n d β 2 = E [ ( X E ( X ) ) 4 ] ( V ( X ) ) 2 = μ 4 4 μ 1 μ 3 + 6 μ 1 2 μ 2 3 μ 1 4 μ 2 μ 1 2 2 ,
where μ 1 , μ 2 , μ 3 and μ 4 were given in Proposition 4. The result follows replacing the corresponding terms.    □
Figure 4 shows the skewness and kurtosis coefficients for the SMMB and SMB distributions in terms of the shape parameter q. Note that both the skewness and kurtosis coefficients become smaller as q increases.

2.3. Order Statistics

Given a random sample of size n of X SMMB ( β , q ) and denoting by X ( j ) the j-th order statistics, j { 1 , , n } . The following proposition presents the pdf for order statistics from the SMMB model.
Proposition 6.
The pdf of X ( j ) is
f X ( j ) ( x ) = 3 n ! ( j 1 ) ! ( n j ) ! x ( 1 β x 2 ) q + 3 2 2 β 3 / 2 x 3 3 B q 2 , 3 2 j 2 F 1 j 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 × 1 2 β 3 / 2 x 3 3 B q 2 , 3 2 2 F 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 n j .
In particular, the pdf of the minimum, X ( 1 ) , is
f X ( 1 ) ( x ) = 2 n β 3 / 2 x 2 ( 1 β x 2 ) q + 3 2 B q 2 , 3 2 1 2 β 3 / 2 x 3 3 B q 2 , 3 2 2 F 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 n 1 ,
and the pdf of the maximum, X ( n ) , is
f X ( n ) ( x ) = 3 n x ( 1 β x 2 ) q + 3 2 2 β 3 / 2 x 3 3 B q 2 , 3 2 n 2 F 1 n 1 q + 3 2 , 3 2 ; 5 2 ; β x 2 .
Proof. 
Since we are dealing with an absolutely continuous model, the pdf of the j-th order statistics is obtained by applying
f X ( j ) ( x ) = n ! ( j 1 ) ! ( n j ) ! f ( x ) F ( x ) j 1 1 F ( x ) n j ,
where F and f are the cdf and pdf of the SMMB distribution.    □

2.4. Entropy

The following results show the Rényi and Shannon entropy of a random variable X SMMB ( β , q ) . Such results are commonly related with a state of disorder or uncertainty, and are used in different fields such as thermodynamics, statistical physics and information theory.
Proposition 7.
The Rényi ( R γ ) entropy for a random variable X S M M B ( β , q ) and values of γ > 0 and γ 1 is as follows
R γ ( X ) = 1 1 γ ( γ 1 ) log ( 2 β ) + ψ 2 γ + 1 2 + ψ γ q + γ 1 2 + γ κ ( q ) ψ ( γ ρ 3 ) ,
where ρ i = q + i 2 , with i = 0 , 3 , ψ ( · ) is the digamma function and κ ( q ) = ψ ( ρ 3 ) ψ ( ρ 0 ) ψ ( 1.5 ) .
Proof. 
Calculating the Rényi entropy by its definition we have:
R γ ( X ) = 1 1 γ log 0 f X ( x ; β , q ) γ d x = 1 1 γ log 0 2 γ β 3 γ / 2 x 2 γ B γ q 2 , 3 2 1 + β x 2 γ ( q + 3 ) 2 d x ,
then considering the change of variable u = β x 2 , we obtain the result.    □
Corollary 3.
Let X S M M B ( β , q ) , the Shannon (S) entropy for random variable X is given by
S ( X ) = ρ 3 ψ ( ρ 3 ) ρ 1 ψ ( ρ 0 ) + ψ ( ρ 0 ) + ψ ( 1.5 ) ψ ( ρ 3 ) ψ ( 1.5 ) log ( 2 β ) ,
where ρ i = q + i 2 , with i = 0 , 1 , 3 .
Proof. 
The Shannon entropy is obtained for the limit case γ 1 in the definition for the Renyi entropy R γ ( X ) . Therefore, applying L’Hopital’s rule, lim γ 1 R γ ( X ) , the result is obtained.    □

3. Inference

In this section, we present the estimators for the parameters β and q of the SMMB distribution obtained by the moments and maximum likelihood methods.

3.1. Moments Estimators

The following proposition presents the moment estimators for the SMMB distribution.
Proposition 8.
Let X 1 , , X n be a random sample from X S M M B ( β , q ) . Then the moments estimators ( β ^ M , q ^ M ) of ( β , q ) for q > 2 are as follows
β ^ M = 2 Γ q ^ M 1 2 π X ¯ Γ q ^ M 2 2 ,
8 X 2 ¯ Γ 2 q ^ M 1 2 3 π X ¯ 2 Γ q ^ M 2 Γ q ^ M 2 2 = 0 ,
where (9) is solved for q ^ M numerically, then the value of q ^ M found is replaced in (8) and β ^ M is obtained.
Proof. 
From Proposition 4 and using the first two equations, and following the moments method, we have
X ¯ = 2 Γ q 1 2 π β Γ q 2 , X 2 ¯ = 3 Γ q 2 2 2 β Γ q 2 .
Solving the first equation above for β we obtain β ^ M given in (8). Substituting β ^ M in the second equation above, we obtain the result given in (9).    □

3.2. Maximum Likelihood Estimator

Let X 1 , . . . , X n be a random sample from X SMMB ( β , q ) . Then the log-likelihood function is
( β , q ) = n log ( 2 β 3 ) + 2 i = 1 n log x i n log B q 2 , 3 2 q + 3 2 i = 1 n log 1 + β x i 2 .
The score equations are given by
β ( q + 3 ) i = 1 n x i 2 1 + β x i 2 = 3 n ,
n ψ q 2 + i = 1 n log ( 1 + β x i 2 ) = n ψ q + 3 2 .
Solutions for Equations (11) and (12) can be obtained by using numerical procedures such as the Newton-Raphson algorithm. An alternative to obtain the ML estimates is by maximizing (10) using the optim subroutine in the R software package [19].

3.3. EM-Algorithm

Employing the stochastic representation of the SMMB model given in Proposition 2, we can apply an iterative method to find maximum likelihood estimators based on the EM algorithm (see Dempster [20]). This will greatly simplify the estimation process, since the steps of the algorithm that we will develop will both have a closed form (both E and M steps). We next present Lemma 2, which will be useful in the application of the EM algorithm.
Lemma 2.
Let X G a m m a ( k , σ ) with k > 0 shape parameter and σ > 0 rate parameter, where density function is f X ( x ) = 1 Γ ( k ) σ k x k 1 e σ x , x > 0 . Then
1 .
X k G G σ m , k m , 1 m , m > 0 with pdf given in (1).
2 .
E [ log ( X ) ] = ψ ( k ) log σ .
In this context, the SMMB distribution can also be written employing the hierarchical representation defined in Proposition 3:
X i Z i = z i M B ( z i 2 β ) a n d Z i G G ( 1 , q , 2 ) , i = 1 , , n .
In this scenario, we have that joint pdf X i and Z i is
f X i , Z i ( x i , z i θ ) = f X i Z i ( x i z i , θ ) · f Z i ( z i θ ) = 8 β 3 / 2 x i 2 z i q + 2 π Γ q 2 exp { z i 2 ( 1 + β x i 2 ) } ,
where θ = ( β , q ) is the vector of parameters. We have that x = ( x 1 , , x n ) and z = ( z 1 , , z n ) represent the observed and unobserved data, respectively. The complete data are given by x c = ( x , z ) . Furthermore, we denote c ( θ x c ) and Q ( θ θ ^ ) = E [ c ( θ | x c ) x , θ ^ ] as the log-likelihood complete function and its expected conditional value over the observed data. From (13), we see directly that the complete log-likelihood function is given by
c ( θ x c ) = i = 1 n log f ( x i , z i | θ ) = ( q + 2 ) i = 1 n log z i + 3 n 2 log β n log Γ q 2 i = 1 n z i 2 ( 1 + β x i 2 ) + C ,
where C is a constant that does not depend on θ . To obtain Q ( θ θ ^ ) , we need to calculate E [ z i 2 x i , θ ] and E [ log ( z i ) x i , θ ] . From Equation (13), it is direct that
f ( z i x i , θ ) z i q + 2 exp { z i 2 ( 1 + β x i 2 ) } z i ( q + 3 ) 1 exp { z i ( 1 + β x i 2 ) 1 / 2 2 } .
Therefore, we conclude that Z i | x i , θ G G ( 1 + β x i 2 ) 1 / 2 , q + 3 , 2 . Then, using Lemma 2 it follows
E [ Z i 2 x i , θ ] = 1 + β x i 2 1 / 2 2 Γ ( q + 3 ) + 2 2 Γ q + 3 2 = q + 3 2 ( 1 + β x i 2 ) , and
E [ log ( Z i ) | x i , θ ] = 1 2 ψ q + 3 2 1 2 log ( 1 + β x i 2 ) .
Thus,
Q ( θ θ ^ ) = ( q + 2 ) i = 1 n log ( z i ) ^ + 3 n 2 log β n log Γ q 2 i = 1 n z i 2 ^ ( 1 + β x i 2 ) ,
where z i 2 ^ = E [ Z i 2 x i , θ = θ ^ ] and log ( z i ) ^ = E [ log ( Z i ) x i , θ = θ ^ ] are given in (15) and (16), respectively. Hence, partially differentiating with respect to β and q and equalising to zero, we obtain the following expressions:
β ^ = 3 n 2 i = 1 n z i 2 ^ x i 2 and q ^ = 2 ψ 1 2 · log ( z ) ^ ¯ ,
where log ( z ) ^ ¯ is the mean of log ( z i ) ^ s and ψ 1 ( · ) is the inverse of the digamma function. It is possible to check that Q ( θ | θ ) ^ / β and Q ( θ | θ ) ^ / q are increasing functions in β and q, respectively. This guarantees that β ^ and q ^ are the unique maximum for this function. Thus, the EM algorithm is reduced to the following steps:
  • Step-E: For i = 1 , , n compute
    z i 2 ^ ( k + 1 ) = q ( k ) + 3 2 1 + β ( k ) x i 2 a n d log ( z i ) ^ ( k + 1 ) = 1 2 ψ q ( k ) + 3 2 log 1 + β ( k ) x i 2 .
  • Step-M: Update the parameters as
    β ^ ( k + 1 ) = 3 n 2 i = 1 n z i 2 ^ ( k + 1 ) x i 2 a n d q ^ ( k + 1 ) = 2 ψ 1 2 · log ( z ) ^ ( k + 1 ) ¯ .
Steps E and M are repeated until we reach a defined convergence criterion. For instance, we specify that the difference between the successively obtained values should be inferior to a pre-established value, i.e.,
max β ( k ) β ( k 1 ) , q ( k ) q ( k 1 ) < ϵ , w i t h ϵ = 10 4 .

3.4. Fisher’s Information Matrix

Let us now consider X S M M B ( β , q ) . For a single observation x of X, the log-likelihood function for θ = ( β , q ) is given by
( θ ) = log ( 2 ) + 3 2 log ( β ) + 2 log ( x ) log B q 2 , 3 2 q + 3 2 log ( 1 + β x 2 ) .
The corresponding first and second partial derivatives of the log-likelihood function are derived in Appendix A. It can be shown that the Fisher’s information matrix, denoted by I F ( · ) , for the S M M B distribution is provided by
I F ( θ ) = 3 q 2 β 2 ( q + 5 ) 3 2 β ( q + 3 ) 3 2 β ( q + 3 ) 1 4 ψ q + 3 2 ψ q 2 ,
where ψ ( · ) is the trigamma function. Hence, for large samples, the ML estimator, θ ^ , of θ is asymptotically normal bivariate, i.e.,
n θ ^ θ D N 2 ( 0 , I F ( θ ) 1 ) , a s   n + .
As a result, the asymptotic covariance matrix of the ML estimators θ ^ is the inverse of Fisher’s information matrix I F ( θ ) . Since the parameters are unknown, the observed information matrix is usually considered, where the unknown parameters are estimated by ML. Note that asymptotic confidence intervals for β and q with confidence level 100 ( 1 α ) % should be constructed as
C I β , 100 ( 1 α ) % = β ^ z α / 2 × s e ^ β ^ a n d C I q , 100 ( 1 α ) % = q ^ z α / 2 × s e ^ q ^ ,
where s e ^ β ^ and s e ^ q ^ are the root of the first and second element of the diagonal of I F ( θ ) 1 , respectively, and z α / 2 satisfies P ( Z > z α / 2 ) = α / 2 , with Z a random variable with standard normal distribution.

4. Simulation Study

In this section, we present a simulation study to assess the performance of the EM algorithm for the estimators of β and q in the SMMB model. We consider 1000 replicates of three sample sizes generated using the SMMB model: n = 50 , 100 and 200. By using the stochastic representation given in Equation (2), it is possible to generate random numbers for the SMMB( β , q ) distribution, which leads to the following algorithm:
1.
Generate P i χ 3 2 (chi squared with 3 degrees of freedom), i = 1 , , n .
2.
Compute V i = P i 2 β , i = 1 , , n .
3.
Generate W i G G ( 1 , q , 2 ) , i = 1 , , n .
4.
Compute X i = V i W i i = 1 , , n .
It then follows that X S M M B ( β , q ) .
For each sample generated, ML estimates were calculated using the EM algorithm. In Table 2 the empirical bias (Bias), the mean of the standard errors (SE), the root of the empirical mean squared error (RMSE) and the 95% coverage probability (CP) based on the asymptotic distribution for ML estimators are given for the parameter estimators. We conclude that the ML estimates have desirable properties. Bias is reasonable and decreases as the sample size increases. As expected, SE and RMSE terms are closer when the sample size increases, which suggests that the SE estimates are well estimated.

5. Application

In this section, we present an illustration using a real dataset in order to compare the fit of the SMMB model with those of the SMB, W and BS distributions. These models are compared using the Akaike Information Criterion (AIC) (see Akaike [21]), and the Bayesian Information Criterion (BIC), introduced in Schwarz [22]. The data set is related to zinc (Zn) concentrations in 86 soil samples obtained from the Mining Department, Universidad de Atacama, Chile; it is available in Reyes et al. [23]. Table 3 provides the descriptive summary for the data, including the sample asymmetry coefficient b 1 and sample kurtosis coefficient b 2 .
Figure 5 and Table 3 show the type of situation in which the SMMB distribution is appropriate: distribution with a very heavy tail (see the typical boxplot), and high sample kurtosis coefficient ( b 2 = 32.3421 ).
Table 4 shows the results for the fit of the BS, W, SMB and SMMB models. Note that the SMMB model provides a better fit than the others since its AIC and BIC values are less than the rest of models. Finally, Figure 6 shows the qq-plot comparing the quantiles for the zinc dataset with the quantiles for the fitted SMMB model, and the empirical cdf for the zinc data compared with the estimated cdf for the fitted SMMB, SMB, BS and W distributions. The first plot suggests that the SMMB model is appropriated for this data and the second confirms that the SMMB model provides a better fit for this data than the other models considered.

6. Conclusions

In this paper, we introduce a new model with a greater kurtosis than the SMB distribution. To estimate the parameters, we use the EM algorithm, obtaining acceptable results for the ML estimators. In the application to a real data set, it shows a better fit than other known models. Some additional characteristics are:
  • SMMB distribution has a more flexible kurtosis coefficient than the SMB distribution, as is clearly shown in Figure 4 (Right panel)
  • Closed expressions are given for its main characteristics: pdf, cdf, moments and coefficients of skewness and kurtosis.
  • We discuss the hazard and survival functions, which are in terms of the hypergeometric function and the order statistics of the SMMB model.
  • Employing the scale mixture representation, the EM algorithm was implemented to calculate the ML estimators.
  • The results of a simulation study indicate that, with a reasonable sample size, an acceptable bias is obtained.
  • An illustration with real data shows that the SMMB model achieves a better fit in terms of the AIC and BIC criteria.
In future research, we plan to extend the proposed model to create a multivariate scale mixture of the Maxwell–Boltzmann model [24] to handle multivariate/clustered and positive data in order to explore other estimation methods for the model [25].

Author Contributions

Conceptualization, J.S.C., K.P.G., H.A.M. and H.W.G.; methodology, H.W.G. and O.V.; software, J.S.C. and D.I.G.; validation, D.I.G., O.V. and M.B.; formal analysis, J.S.C., K.P.G., H.A.M. and D.I.G.; investigation, M.B. and H.W.G.; resources, H.W.G. and O.V.; data curation, J.S.C.; writing—original draft preparation, J.S.C., K.P.G. and H.A.M.; writing—review and editing, D.I.G., M.B., O.V. and H.W.G.; supervision, D.I.G. and H.W.G. All authors have read and agreed to the published version of the manuscript.

Funding

The research of J.S.C., K.P.G., H.A.M. and H.W.G. ómez was supported by Semillero UA-2022.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data available in Reyes et al. [23].

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The first derivatives of ( θ ) are given by
( θ ) β = 3 2 β q + 3 2 x 2 1 + β x 2 , ( θ ) q = 1 2 ψ q 2 1 2 ψ q + 3 2 1 2 log ( 1 + β x 2 ) .
The second derivatives of l ( θ ) are:
2 ( θ ) β 2 = 3 2 β 2 + q + 3 2 x 4 ( 1 + β x 2 ) 2 , 2 ( θ ) β q = x 2 2 ( 1 + β x 2 ) , 2 ( θ ) q 2 = 1 4 ψ q 2 1 4 ψ q + 3 2 ,
where ψ ( · ) and ψ ( · ) are the digamma and trigamma functions respectively.

References

  1. Maxwell, J. On the dynamical theory of gases. Philos. Trans. R. Soc. Lond. 1867, 35, 185–217. [Google Scholar]
  2. Tyagi, R.K.; Bhattacharya, S.K. Bayes estimation of the Maxwell’s velocity distribution function. J. Stat. Comput. Simul. 1989, 29, 563–567. [Google Scholar]
  3. Bekker, A.; Roux, J. Reliability characteristics of the Maxwell distribution: A bayes estimation study. Commun. Stat. Theory Methods 2005, 34, 2169–2178. [Google Scholar] [CrossRef]
  4. Sharma, V.K.; Bakouch, H.S.; Suthar, K. An extended Maxwell distribution: Properties and applications. Commun. Stat. Simul. Comput. 2017, 46, 6982–7007. [Google Scholar] [CrossRef]
  5. Vivekanand, H.K.; Kumar, K. Estimation in Maxwell distribution with randomly censored data. J. Stat. Comput. Simul. 2015, 85, 4264–4274. [Google Scholar] [CrossRef]
  6. Iriarte, Y.A.; Astorga, J.M.; Gómez, H.W. Gamma-Maxwell distribution. Commun. Stat. Theory Methods 2017, 46, 4264–4274. [Google Scholar] [CrossRef]
  7. Dey, S.; Dey, T.; Ali, S.; Mulekar, M.S. Two-parameter Maxwell distribution: Properties and different methods of estimation. J. Stat. Theory. Pract. 2016, 10, 394–403. [Google Scholar] [CrossRef]
  8. Sharma, V.K.; Dey, S.; Singh, S.K.; Manzoor, U. On length and area-biased Maxwell distributions. Commun. Stat. Simul. Comput. 2017, 47, 1506–1528. [Google Scholar] [CrossRef]
  9. Segovia, F.; Gómez, Y.M.; Gallardo, D.I. Exponentiated power Maxwell distribution with quantile regression and applications. SORT 2021, 45, 181–200. [Google Scholar]
  10. Shakil, M.; Golam, B.; Chang, K. Distributions of the product and ratio of Maxwell and Rayleigh random variables. Stat. Pap. 2008, 49, 729–747. [Google Scholar] [CrossRef]
  11. Acitas, S.; Arslan, T.; Senoglu, B. Slash Maxwell distribution: Definition, modified maximum likelihood estimation and applications. Gazi Univ. J. Sci. 2020, 33, 249–263. [Google Scholar] [CrossRef]
  12. Reyes, J.; Barranco-Chamorro, I.; Gómez, H.W. Generalized modified slash distribution with applications. Commun. Stat.-Theory Methods 2020, 49, 2025–2048. [Google Scholar] [CrossRef]
  13. Astorga, J.M.; Reyes, J.; Santoro, K.I.; Venegas, O.; Gómez, H.W. A Reliability Model Based on the Incomplete Generalized Integro-Exponential Function. Mathematics 2020, 8, 1537. [Google Scholar] [CrossRef]
  14. Stacy, E.W. A generalization of the gamma distribution. An. Math. Stat. 1962, 33, 1187–1192. [Google Scholar] [CrossRef]
  15. Abramowitz, M.; Stegun, I.E. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th ed.; National Bureau of Standars: Washington, DC, USA, 1970. [Google Scholar]
  16. Singh, P. Hypergeometric Functions on Cumulative Distribution Function. Asian Res. J. Math. 2019, 13, 1–11. [Google Scholar] [CrossRef] [Green Version]
  17. Glaser, R.E. Bathtub and Related Failure Rate Characterizations. J. Am. Stat. Assoc. 1980, 75, 667–672. [Google Scholar] [CrossRef]
  18. Weisstein, E.W. “Hypergeometric Function”. From MathWorld—A Wolfram Web Resource. Available online: https://mathworld.wolfram.com/HypergeometicFunction.html (accessed on 2 January 2023).
  19. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2022; Available online: https://www.R-project.org/ (accessed on 2 January 2023).
  20. Dempster, A.P.; Laird, N.M.; Rubim, D.B. Maximum likelihood from incomplete data via the EM algorithm (with discussion). J. R. Stat. Soc. Ser. B 1977, 39, 1–38. [Google Scholar]
  21. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control. 1974, 19, 716–723. [Google Scholar] [CrossRef]
  22. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978, 6, 461–464. [Google Scholar] [CrossRef]
  23. Reyes, J.; Rojas, M.A.; Venegas, O.; Gómez, H.W. Nakagami Distribution with Heavy Tails and Applications to Mining Engineering Data. J. Stat. Theory Pract. 2020, 14, 55. [Google Scholar] [CrossRef]
  24. Morán-Vásquez, R.A.; Mazo-Lopera, M.A.; Ferrari, S.L.P. Quantile modeling through multivariate log-normal/independent linear regression models with application to newborn data. Biom. J. 2021, 63, 1290–1308. [Google Scholar] [CrossRef] [PubMed]
  25. Lange, K.; Sinsheimer, J.S. Normal/Independent Distributions and Their Applications in Robust Regression. J. Comput. Graph. Stat. 1993, 2, 175–198. [Google Scholar]
Figure 1. (Left panel): pdf for SMMB ( β = 1 , q ) model. (Right panel): cdf for SMMB ( β = 1 , q ) model. In both cases, different values for q were considered.
Figure 1. (Left panel): pdf for SMMB ( β = 1 , q ) model. (Right panel): cdf for SMMB ( β = 1 , q ) model. In both cases, different values for q were considered.
Mathematics 11 00529 g001
Figure 2. (Left panel): survival function for SMMB ( β = 1 , q ) model. (Right panel): hazard function for SMMB ( β = 1 , q ) model. In both cases, different values for q were considered.
Figure 2. (Left panel): survival function for SMMB ( β = 1 , q ) model. (Right panel): hazard function for SMMB ( β = 1 , q ) model. In both cases, different values for q were considered.
Mathematics 11 00529 g002
Figure 3. f q ( x ) function graphs for different values of q.
Figure 3. f q ( x ) function graphs for different values of q.
Mathematics 11 00529 g003
Figure 4. (Left panel): skewness coefficient. (Right panel): kurtosis coefficient; for the SMMB and SMB distributions.
Figure 4. (Left panel): skewness coefficient. (Right panel): kurtosis coefficient; for the SMMB and SMB distributions.
Mathematics 11 00529 g004
Figure 5. Boxplot for zinc dataset.
Figure 5. Boxplot for zinc dataset.
Mathematics 11 00529 g005
Figure 6. (Left panel): qq-plot for zinc data set versus the fitted SMMB distribution. (Right panel): Empirical cdf (black) and estimated cdf of SMMB (red), SMB (blue), BS (green) and W (orange) models for the zinc data set.
Figure 6. (Left panel): qq-plot for zinc data set versus the fitted SMMB distribution. (Right panel): Empirical cdf (black) and estimated cdf of SMMB (red), SMB (blue), BS (green) and W (orange) models for the zinc data set.
Mathematics 11 00529 g006
Table 1. Specific values of 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 ) for different q values.
Table 1. Specific values of 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 ) for different q values.
q 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 ) q 2 F 1 ( 3 2 , q + 3 2 ; 5 2 ; β t 2 )
1 3 2 β 3 / 2 t 3 tan 1 ( β t ) β t 1 + β t 2 5 1 16 β 3 / 2 t 3 3 tan 1 ( β t ) + β t ( β t 2 + 3 ) ( 3 β t 2 1 ) ( 1 + β t 2 ) 3
2 1 / ( 1 + β t 2 ) 3 / 2 6 4 β t 2 ( 2 β t 2 + 7 ) + 35 35 ( 1 + β t 2 ) 7 / 2
3 3 8 β 3 / 2 t 3 tan 1 ( β t ) + β t ( β t 2 1 ) ( 1 + β t 2 ) 2 7 1 128 β 3 / 2 t 3 β t ( β t 2 ( 5 β t 2 ( 3 β t 2 11 ) + 73 ) 15 ) ( 1 + β t 2 ) 4 + 15 tan 1 ( β t )
4 5 + 2 β t 2 5 ( 1 + β t 2 ) 5 / 2 8 105 + 2 β t 2 ( 4 β t 2 ( 2 β t 2 + 9 ) + 63 ) 105 ( 1 + β t 2 ) 9 / 2
Table 2. Simulation study for different combinations of parameters for the SMMB( β , q ) model.
Table 2. Simulation study for different combinations of parameters for the SMMB( β , q ) model.
True Value n = 50 n = 100 n = 200
β q Estim.BiasSERMSECPBiasSERMSECPBiasSERMSECP
31 β ^ 0.1131.2941.33490.10.0740.8970.92592.10.0270.6230.66192.2
q ^ 0.0570.2160.23296.50.0240.1450.15395.00.0160.1020.11194.1
2 β ^ 0.0651.2901.33590.10.0390.8960.92192.50.0190.6280.66591.8
q ^ 0.2110.6130.70995.90.0900.3880.42695.50.0520.2650.29694.9
3 β ^ 0.0401.3961.44489.60.0150.9660.99692.10.0090.6790.70992.0
q ^ 0.5241.3291.62595.70.2260.7580.89295.10.1150.4950.55895.2
51 β ^ 0.1882.1552.22290.10.1221.4951.54192.10.0441.0381.10292.2
q ^ 0.0570.2160.23296.50.0240.1450.15395.10.0160.1020.11194.1
2 β ^ 0.1062.1492.22490.10.0641.4931.53492.50.0301.0461.10891.8
q ^ 0.2110.6130.70995.90.0910.3880.42695.50.0530.2650.29694.9
3 β ^ 0.0612.3242.41089.50.0241.6101.65992.10.0141.1321.18292.0
q ^ 0.5671.4782.12095.70.2270.7580.89295.10.1160.4950.55895.2
71 β ^ 0.2623.0173.11190.10.1692.0932.15792.10.0611.4531.54292.2
q ^ 0.0570.2160.23296.50.0240.1450.15395.10.0160.1020.11194.1
2 β ^ 0.1483.0083.11390.10.0882.0902.14892.50.0411.4641.55091.8
q ^ 0.2110.6130.70995.90.0910.3880.42695.50.0530.2650.29694.9
3 β ^ 0.0783.2513.37889.40.0332.2542.32292.10.0191.5851.65592.0
q ^ 0.6101.5662.51995.70.2270.7580.89295.10.1160.4950.55995.2
Table 3. Descriptive statistics for zinc data set.
Table 3. Descriptive statistics for zinc data set.
n x ¯ s b 1 b 2
8696.721148.4345.08832.342
Table 4. Estimates, SE in parenthesis, log-likelihood, AIC and BIC values for zinc concentration data.
Table 4. Estimates, SE in parenthesis, log-likelihood, AIC and BIC values for zinc concentration data.
ParametersBSWSMBSMMB
α 1.3038 (0.0995)0.0125 (0.0047)3.4666 (0.1226)-
β 50.8841 (5.8694)0.9632 (0.0701)-0.0007 (0.0002)
q--0.4077 (0.1827)1.6506 (0.3079)
log-likelihood−484.7785−479.0418−471.7718−469.5580
AIC973.557962.084947.544943.1161
BIC978.466966.992952.452948.0248
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Castillo, J.S.; Gaete, K.P.; Muñoz, H.A.; Gallardo, D.I.; Bourguignon, M.; Venegas, O.; Gómez, H.W. Scale Mixture of Maxwell-Boltzmann Distribution. Mathematics 2023, 11, 529. https://doi.org/10.3390/math11030529

AMA Style

Castillo JS, Gaete KP, Muñoz HA, Gallardo DI, Bourguignon M, Venegas O, Gómez HW. Scale Mixture of Maxwell-Boltzmann Distribution. Mathematics. 2023; 11(3):529. https://doi.org/10.3390/math11030529

Chicago/Turabian Style

Castillo, Jaime S., Katherine P. Gaete, Héctor A. Muñoz, Diego I. Gallardo, Marcelo Bourguignon, Osvaldo Venegas, and Héctor W. Gómez. 2023. "Scale Mixture of Maxwell-Boltzmann Distribution" Mathematics 11, no. 3: 529. https://doi.org/10.3390/math11030529

APA Style

Castillo, J. S., Gaete, K. P., Muñoz, H. A., Gallardo, D. I., Bourguignon, M., Venegas, O., & Gómez, H. W. (2023). Scale Mixture of Maxwell-Boltzmann Distribution. Mathematics, 11(3), 529. https://doi.org/10.3390/math11030529

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop