Parameter Estimation of the Extended Generalized Gaussian Family Distributions using Maximum Likelihood Scheme

INTRODUCTION

In the field of science and engineering, the family of the Generalized Gaussian Distributions (GGD) are applied to model and detect the Gaussian and non-Gaussian signals (Miller and Thomas, 1972; Kassam, 1988). For example, in natural environments, the performance of communication systems can significantly degrade due to the atmospheric noise caused by the lightning activity in local and distant storm, or the ambient noise in the ocean. In seismic research (Walden, 1985), the non-Gaussian distribution has been used to model the amplitude distribution of the primary reflection coefficients. The aim is to estimate the reflection coefficient sequences as accurately as possible. In image processing, the fundamental work is fitting the image intensity histogram data to gain underlying knowledge of data source. Fan and Lin (2009) developed an expectation maximization algorithm to estimate the parameter of the GGD mixture model that best fit the observed image intensity histogram data. Other useful applications of the GGD also include, video coding (Sharifi and Leon-Garcia, 1995), handwritten character identification (Solihin and Leedham, 1999), speech recognition (Gazor and Zhang, 2003), etc. Therefore, the choice between different families of distribution is important and the estimation of the designed probability density function from observed data is the major task.

In this study, the parameter estimator for the extended generalized Gaussian distribution to include asymmetry is applied. The parameter estimation is determined by the maximum likelihood scheme with iteration method. The key factor of the iteration method depends on the choice of the initial and step values. Soderlind (2006) presented and discussed the approaches to step size selection in control theory and signal processing. Tlelo-Cuautle et al. (2007) proposed a procedure to control the step-size, where the step size is determined by the eigenvalues of its state matrixes. For saving the computing cost and improving the accuracy, we proposed a table look-up method.

THE MODEL

The generalized Gaussian distribution to include asymmetry proposed by Nandi and Mampel (1995) is:

(1)

and

(2)

where, α_n>0, α_p>0, β_n>0 and Γ is gamma function. By setting:

(3)

then the variance of Eq. 1 is chosen as one. By changing the parameters in Eq. 1, one can obtain symmetric (α_n = α_p) and asymmetric (α_n ≠ α_p) distributions. For symmetric cases, Fig. 1, (α_n , α_p) are (1,1) and (2,2) represent the Laplace distribution and Gaussian distribution, respectively. When α_n tends to 0 the distribution tends to impulse and when α_n tends to infinite the distribution tends to Uniform distribution. For asymmetric cases, (Fig. 2), the right skew distribution can be modeled by reducing the parameter α_p and as α_p decreases the tail length in the right side increases. Therefore, Eq. 1 provides a wide range of uni-modal distributions with symmetry and asymmetry.

Fig. 1:	Symmetric model for varying the parameters of (αn = αp)

Fig. 2:	Asymmetric model for varying the parameters of (αn = αp)

THE ESTIMATOR

Given a series samples x_i, where i = 1,2,3,...,N, the likelihood function of Eq. 1 is:

(4)

where, N₁ (N_r) is the number of x_i samples <0 (≥0). Taking the logarithm of Eq. 4, one gets:

(5)

If ∂ InL/∂α_n, ∂InL/∂β_n and ∂InL/∂α_p are equate to zero, then the maximum likelihood estimates α_n, β_n and α_p are obtained, i.e.,

(6)

(7)

(8)

Where,

(9)

and Ψ(•) is digamma function defined as:

(10)

To fit the observed data and approximate the solutions α_n, β_n and α_n, one need to solve (Eq. 6-8). In the following two iterated methods are developed.

Iterated method 1: Given a reasonable initial values α_n0, β_n and α_p, (Eq. 6-8) can be solved iteratively for α_n, β_n and α_p. First step, compute Eq. 6 with (α_n0±Δα_n) by fixed β_n0 and α_p0, then α_n1 is selected by the value of (α_n0±Δα_n) or (α_n0-Δα_n) so that the likelihood function of Eq. 6 approaches zero. Second step, compute Eq. 7 with (β_n0±Δβ_n) by fixed α_n1 and α_p0, then β_n1 can be updated by the value of (β_n0±Δβ_n) or (β_n0-Δβ_n) so that Eq. 7 approaches zero. Third step, compute Eq. 8 with (α_n0±Δα_n) by fixed α_n1 and β_n1, then α_p1 can be updated by the value of (α_n0+Δα_n) or (α_n0-Δα_n) so that Eq. 8 approaches zero. Thus, one iteration is completed for updating α_n1, β_n1 and α_n1. This process continues until Eq. 6-8 simultaneously approach zero with some error tolerances. The accuracy of the estimation of this method depends on the choice of the step values Δα_n, Δβ_n and Δα_p. Therefore, small step value leads to accuracy of the estimation but requires more computing cost.

Iterated method 2: In this section the gradient method for iteration is developed. In the j-th iteration the Taylor series of Eq. 11 is applied to update the parameters α_nj, β_nj and α_pj.

(11)

where, θ = (α_n, β_n, α_p)^T, grad In L = (∂In L/∂α_n, ∂In L/∂β_n, ∂In L/∂α_p)^T and:

(12)

is a Fisher information matrix. The inverse of the Fisher information matrix plays a key role for the iteration of Eq. 11. Given N samples of the observed data, the Fisher matrix can be obtained from Eq. 12. For each iteration, the expected value of the information matrix is calculated with updating parameters α_nj, β_nj and α_pj. The estimation of the parameter is obtained as the value V(θ^(j))^-1 converge to zero with some error tolerances.

Cramer-Rao lower bound: It is important to understand the unbiasedness of an estimator. The most common way to decide whether or not an estimator’s distribution is unbiased is by looking at its expected value. The Cramer-Rao Lower Bound (CRLB) for the standard deviation of every unbiased estimator is defined as (Hogg and Tanis, 1989):

(13)

(14)

(15)

where, N is the samples of the every realizations and the E[•] is the expected value. Note that Eq. 13-15 are equivalence to the inverse of the square root of the diagonal elements in Eq. 12.