# Mixture of Generalized Gamma Density-Based Score Function for Fastica

Authors: M. EL-Sayed Waheed; Osama Abdo Mohamed; M. E. Abd El-aziz
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We propose an entirely novel family of score functions for blind signal separation (BSS), based on the family of mixture generalized gamma density which includes generalized gamma, Weilbull, gamma, and Laplace and Gaussian probability density functions. To blindly extract the independent source signals, we resort to the FastICA approach, whilst to adaptively estimate the parameters of such score functions, we use Nelder-Mead for optimizing the maximum likelihood (ML) objective function without relaying on any derivative information. Our experimental results with source employing a wide range of statistics distribution show that Nelder-Mead technique produce a good estimation for the parameters of score functions.


---

## Body

## 1. Introduction

By definition, independent component analysis (ICA) is the statistical method that searches for a linear transformation, which can effectively minimize the statistical dependence between its components [1]. Under the physically plausible assumption of mutual statistical independence between these components, the most application of ICA is blind signal separation (BSS). In its simplest form, BSS aims to recover a set of unknown signals, the so-called original sources  s(t)=[s1(t),s2(t),…,sn(t)]T∈Rn, by relying exclusively on information that can be extracted from their linear and instantaneous mixtures x(t)=[x1(t),x2(t),…,xm(t)]T∈Rm, given by (1.1)x(t)=As(t),t=1,2,…,m,
where A∈Rmxn is an unknown mixing matrix of full rank and m≥n. In doing so, BSS remains truly (blind)   in the sense that very little to almost nothing be known a priori for the mixing matrix or the original source signals.Often sources are assumed to be zero-mean and unit-variance signals with at most one having a Gaussian distribution. The problem of source estimation then boils down to determining the unmixing matrixW∈Rnxm such that the linear transformation of the sensor observation is (1.2)u(t)=Wx(t),t=1,2,…,n,
where u(t)=[u1(t),u2(t),…,un(t)]T∈Rn yield an estimate of vector s(t) corresponding to the original or true sources. In general, the majority of BSS approaches perform ICA, by essentially optimizing the negative log-likelihood (objective) function  with  respect  to  the  unmixing  matrix W such that  (1.3)L(u,W)=∑l=1nE[log⁡pul(ul)]-log⁡|det⁡(W)|,
where E[·] represents the expectation operator and pul(ul) is the model for the marginal probability density function (pdf) of ul, for all l=1,2,…,n. Normally, matrix W is regarded as the parameter of interest and the pdfs of the sources are considered to be nuisance parameters. In effect, when correctly hypothesizing upon the distribution of the sources, the maximum likelihood (ML) principle leads to estimating functions, which in fact are the score functions of the sources [2](1.4)φl(ul)=-ddullog⁡pul(ul).
In principle, the separation criterion in (1.3) can be optimized by any suitable ICA algorithm where contrasts are utilized (see; e.g., [2]). A popular choice of such a contrast-based algorithm is the so-called fast (cubic) converging Newton-type (fixed-point) algorithm, normally referred to as FastICA [3], based on(1.5)Wk+1=Wk+D(E[φ(u)uT]-diag⁡(E[φl(ul)ul]))Wk,
where, as defined in [4],(1.6)D=diag⁡(1(E[φl(ul)ul]-E[φl′(ul)])),
with φ(t)=[φ1(u1),φ2(u2),…,φn(un)]T being valid for all l=1,2,…,n. In the ICA framework, accurately estimating the statistical model of the sources at hand is still an open and challenging problem [2]. Practical BSS scenarios employ difficult source distributions and even situations where many sources with very different pdfs are mixed together. Towards this direction, a large number of parametric density models have been made available in recent literature. Examples of  such  models  include the generalized Gaussian density (GGD) [5], the generalized lambda density (GLD), and the generalized beta distribution (GBD) or even combinations and generalizations such as super and generalized Gaussian mixture model (GMM) [6], the generalized gamma density (GGD) [7], the Pearson family of distributions [4], and even the so-called extended generalized lambda distribution (EGLD) which is an extended parameterizations of the aforementioned GLD and GBD models [8]. In the following section, we propose Mixture Generalized Gamma Density (MGΓD) for signal modeling in blind signal separation.


## 2. Mixture Generalized Gamma Density (MGΓD)

A Mixture Generalized Gamma Density (MGΓD) is a parametric statistical model which assumes that the data originates from weighted sum of several generalized gamma sources [9]. More specifically, a MGΓD is defined as(2.1)p(u∣θ)=∑i=1kmipi(u∣ci,αi,βi,γi),
where(i)
θ=(mi,ci,αi,βi,γi),i=1,2,…,k,(ii)
k  is the number of mixture density components,(iii)
mi is the ith mixture weight and satisfies mi≥0,∑i=1kmi=1,(iv)
pi(u|ci,αi,βi,γi) is an individual density of the generalized gamma density which is characterized by [7],(2.2)pi(u∣ci,αi,βi,γi)=γiβiαi|u-ci|αiγi-12Γ(αi)exp⁡(-βi|u-ci|γi),
where ci is the location parameter, βi>0 is the scale parameter, αi>0 is the shape/power parameter, and γi>0 is the shape parameter. Γ(α) is the gamma function, defined by(2.3)Γ(x)=∫0∞tx-1e-tdt.
By varying the parameters, it is possible to characterize a large class of distributions such as Gaussian, sub-Gaussian (more peaked, than Gaussian, heavier tail), and supergaussian (flatter, more uniform). It is noticed that for γ=1, the (GΓD) define gamma density as special case. Furthermore, if  γ=2 and α=0.5, it become the Gaussian pdf, and if γ=1and α=1, it represent the Laplacian pdf.Figures1 and 2 show some examples of pdf for MGΓD for k=1 and k=2. Thanks to the shape parameters, the MGΓD is more flexible and can approximate a large class of statistical distributions, this distribution requires to estimate 5*k parameters, θ(mi,ci,αi,βi,γi)i=1,2,…,k. Particularly, we discuss the estimation of these parameters in detail in the following section.Figure 1
The probability density function of MGΓD fork=1, c=3, α=0.5, β=3, and  γ=6,2,3.5,1(p1,p2,p3,p4) respectively. The MGΓD allows us to include the case of Gaussian mixture model α=0.5 and γ=2.Figure 2
The probability density function of MGΓD fork=2,  c1=2, c2=7,  α1=α2=0.5, and β1=β2=2 with different weight m1=0.7,0.3,0.4,0.6, m2=1-m1, and power parameters γ1=1.9,2,2.5,3 and γ2=2,4,4.3,1.5(g,f,h,w), respectively.


## 3. Numerical Optimization of the Log-Likelihood Function to Estimate MGΓD Parameters

We propose in this section a generalization of the method proposed in [9] which address only the case of 2 components by setting the derivatives of the log-likelihood function to zeros. The log-likelihood function of (2.2), given by [10](3.1)L(u∣θ)=∑i=1k∑j=1N-1hi,jlog⁡[mipi(uj∣ci,αi,βi,γi)],
where N the sample size and hi,j=p(i∣uj)(i=1,…,k,j∈[0,N-1]) represents the conditional expectation of pi given the observation uj, this means the posterior probability that uj belongs to the ith component. In the case of generalized gamma distribution, if we substitute (2.2) into (3.1) and after some manipulation, we obtain the following form of L(u∣θ)(3.2)L(U∣θ)=∑i=1k∑j=1Nhi,jlog⁡(mi)+∑i=1k∑j=1Nhi,j[log⁡(γi)+(αiγi-1)log⁡|uj-ci|-log⁡(2)+αilog⁡(βi)-log⁡(Γ(αi))-βi|uj-ci|γi].
Accordingly, we obtain for i=1,2,…,k the following nonlinear equation related to the estimated parameters by derivatives of the log-likelihood function with respect to ci,αi,βi, and  γiand setting these derivatives to zeros, we obtain (3.3)∂L(U∣θ)∂ci=0,∂L(U∣θ)∂βi=0,∂L(U∣θ)∂γi=0,∂L(U∣θ)∂αi=0,ψ(αi)=log⁡(βi)+∑j=1Nhijlog⁡|uj-ci|γi∑j=1N-1hij,βi=αi∑j=1Nhij∑j=1Nhij|uj-ci|i,(3.4)∑j=1Nhij+αi∑j=1Nhijlog⁡|uj-ci|γi=βi∑j=1Nhij|uj-ci|γilog⁡|uj-ci|γi,
where Ψ(·) is the digamma function (Ψ(x)=Γ′(x)/Γ(x)). After a little mathematical manipulation, the ML estimate of γi is obtained(3.5)ψ((∑j=1Nhij)(∑j=1Nhij|uj-ci|γi)(∑j=1Nhij)(∑j=1Nhij|uj-ci|γilog⁡|uj-ci|γi)-(∑j=1Nhij|uj-ci|γi∑k=1Nhiklog⁡|uk-ci|γi))-log⁡((∑j=1Nhij)2(∑j=1Nhij)(∑j=1Nhij|uj-ci|γilog⁡|uj-ci|γi)-(∑j=1Nhij|uj-ci|γi∑k=1Nhiklog⁡|uk-ci|γi))-∑j=1Nhijlog⁡|uj-ci|γi∑j=1Nhij=0.
Given the estimate of γi, it is straightforward to derive the estimate for αi,βi, and ci. Let γ̂i be the estimate of γi. Then, (3.6)α̂i=(∑j=1Nhij)(∑j=1Nhij|uj-ci|γi)(∑j=1Nhij)(∑j=1Nhij|uj-ci|γilog⁡|uj-ci|γi)-(∑j=1Nhij|uj-ci|γi∑k=1Nhiklog⁡|uk-ci|γi),β̂i=α̂i∑j=1Nhij∑j=1Nhij|uj-ci|γi,∑j=1Nhi,j[(α̂iγ̂i-1)η(uj)-γ̂iβ̂iη(uj)|uj-ĉi|γ̂i]=0,
where (3.7)η(x)={1ifuj-ci<0,-1ifuj-ci≥0,
where α̂i and β̂i are the resulting estimates for αi and βi, respectively, and to estimate the location parameter, we solve (3.4) by gradient ascent. The estimation of weight coefficient obtains directly from hi,j as follows [10](3.8)m̂i=1N∑j=1Nhi,j.However, (3.5) cannot be easily solved, so we adopt the gradient ascent algorithm to obtain the estimate of γi and determine the estimates of αi, βi, and ci uniquely using this value of γi.Alternative numerical method can be used to estimate the parameters is called NM, where the appeal of the NM optimization technique lies in the fact that it can minimize the negative of the log-likelihood objective function given in (3.2), essentially without relying on any derivative information. Despite the danger of unreliable performance (especially in high dimensions), numerical experiments have shown that the NM method can converge to an acceptably accurate solution with substantially fewer function evaluations than multidirectional search descent methods. Good numerical performance and a significant improvement in computational complexity for our estimation method, therefore, optimization with the NM technique, produce a good estimation for parameters in MGΓD. To show the performance of NM, we consider the next example.


### 3.1. Example

We generate random number from MGΓD with parametersk=2, m1=0.25, m2=0.75, α1=0.5, α2=0.5,  β1=2, β2=2, γ1=1, γ2=4, c1=0, and c2=10. By performs NM, we obtain best estimation for parameters. As we show in Table 1, the first 5th values of estimated parameters after being sorted according the value of function. In the following section, we resolve to FastICA algorithm for blind signal separation (BSS), this algorithm depends on the estimated parameters and an unmixing matrix W which estimated by FastICA algorithm.Table 1Estimatedparametersm1c1α1β1γ1m2c2α2β2γ21st  value0.2400.481.860.950.769.740.461.863.962nd value0.2400.471.700.870.769.680.491.903.893rd value0.2300.481.990.900.769.700.491.893.704th value0.2300.471.800.940.779.800.491.703.805th  value0.2200.491.760.950.789.840.471.833.98


## 3.1. Example

We generate random number from MGΓD with parametersk=2, m1=0.25, m2=0.75, α1=0.5, α2=0.5,  β1=2, β2=2, γ1=1, γ2=4, c1=0, and c2=10. By performs NM, we obtain best estimation for parameters. As we show in Table 1, the first 5th values of estimated parameters after being sorted according the value of function. In the following section, we resolve to FastICA algorithm for blind signal separation (BSS), this algorithm depends on the estimated parameters and an unmixing matrix W which estimated by FastICA algorithm.Table 1Estimatedparametersm1c1α1β1γ1m2c2α2β2γ21st  value0.2400.481.860.950.769.740.461.863.962nd value0.2400.471.700.870.769.680.491.903.893rd value0.2300.481.990.900.769.700.491.893.704th value0.2300.471.800.940.779.800.491.703.805th  value0.2200.491.760.950.789.840.471.833.98


## 4. Application of MGΓD in Blind Signal Separation

Novel flexible score function is obtained, by substituting (2.1) into (1.4) for the source estimates ul,l=1,2,…,n, it quickly become obvious that our proposed score function inherits a generalized parametric structure, which in turn can be attributed to the highly flexible MGΓD parent model. In this case, a simple calculus yield the flexible BSS score function (4.1)φl(ul∣θ)=(-1pul(ul))∑i=1kmiγiβiαisign⁡(ul-ci)|ul-ci|(αiγi-2)exp⁡-(βi|ul-ci|γi)2Γ(αi)×[(αiγi-1)-βiγi|ul-ci|γi].
In principle φl(ul∣θ) is capable of modeling a large number of signals, such as speech or communication signals, as well as various other types of challenging heavy- and light-tailed distributions.This is due to the fact that its characterization depends explicitly on all parametersmi,ci,αi,βi,γi, i=1,2,…,k. Other commonly used score functions can be obtained by substituting appropriate values for mi,ci,αi,βi, and γi in (4.1), for instance, when k=1, we have score function (4.2)φl(ul∣θ)=sign⁡(ul-c1)|ul-c1|(β1γ1|ul-c1|γ1-α1γ1+1).
When α1γ1=1 and β1=1, we have a scaled form of the GGD-based score function constitutes such a special case of(4.2)(4.3)φl(ul∣θ)=γ1sign⁡(ul-c1)|ul-c1|γ1-1.
The function φl(ul∣θ) could become singular, in some special cases, essentially those corresponding to heavy-tailed (or sparse) distribution defined for αiγi≥1 with αi=1 and |ul-ci|=0. In practice, to circumvent such deficiency, the denominator in (4.1) can be modified slightly to read
(4.4)φl(ul∣θ)=(-1pul(ul)+ε)∑i=1kmiγisign⁡(ul-ci)|ul-ci|(αiγi-2)exp⁡-(|ul-ci|/βi)γi2βiαiγiΓ(αi)×[(αiγi-1)-γiβiγi|ul-ci|γi],
where ε is a small positive parameter (typically around 104) which, when put to use, can almost always guarantee that the discontinuity of (4.1) or values in or approaching the region |ul-ci|=0 is completely avoided.


## 5. Numerical Experiments

To investigate the separation performance of the proposed MGΓD-based FastICA BSS method, a set of numerical experiments are preformed, in which we consider only two cases  whenk=1, k=2, and we  illustrate this in the following two examples.


### 5.1. Example  1

In this example,k=1 and the data set used consists of different realizations of independent signals, with distributions shown in Table 2. Note that this is a large-scale and substantially difficult separation problem, since it involves a Gaussian, various super- and sub-Gaussian symmetric PDFs, as well as asymmetric distributions. In all cases, the number of data samples has also been designed to be relatively small; for example, N=250. The source signals are mixed (noiselessly) with randomly generated full-rank mixing matrices A. The FastICA method is implemented in the so-called simultaneous separation mode whereas the stopping criterion is set to ε=10-4. FastICA is executed using the flexible MGΓD model is used to model the distribution of the unknown sources, while (3.3)–(3.5) are employed to adaptively calculate the necessary parameters of the MGΓD-based score function defined in (4.4).Table 2SourcesMGΓD(α,γ)PDF plotGaussian0.5, 2.0Laplace1.0, 1.0Uniform0.1, 10Gamma0.5, 1.0Now, to show the performance in this case, we consider three source signals (source), where these signals are generated randomly from Weilbull, gamma, and exponential distributions as follows:(5.1)s1=wblrnd(1:250,1),s2=exprnd(1,1,250),s3=gamrnd(1,1,250).
Let the mixing matrix A and unmixing matrix W be defined as follows:(5.2)A=(.56.79.37.75.65.86.17.32.48),W=(0.2128-0.21770.25660.2184-0.19240.16461.9485-1.87981.8687).
By using the equation x=As, we obtain mixed signals as show in Figure 3, where mixing signals are in the left and source signals are in the right.Figure 3After using FastICA, we recover the sources, and we show the estimated signals in the left and original signals in the right in Figure4 with different in scale only.Figure 4


### 5.2. Example  2

In this example,k=2 in which the data set used consists of different realizations of independent signals. Note that each signal not only a Gaussian, super-, and sub-Gaussian PDFs, but it is mixed of this PDFs as shown in Table 3.Table 3SourcesMGΓD(α1,γ1,α2,γ2)PDF plotGamma_Lapalce0.5, 1.0, 1.0, 1.0Gaussian_Laplace0.5, 2.0, 1.0, 1.0Weilbull_Gamma2.0, 4.0, 0.5, 1.0Uniform_Laplace0.1, 10, 1.0, 1.0To show the performance in this case, we consider three source signals (source) where these signals are generated randomly from Gamma_Lapalce (s1), Weilbull (s2), and Gaussian_Laplace (s3)  distributions. Let the mixing matrix A and unmixing matrix W be defined as follows:(5.3)A=(.56.79.37.75.65.86.17.32.48),W=(-0.12932.6579-4.2322-1.4418-3.15339.9930-2.33108.4107-11.6627).
By using the equation X=AS,we obtain mixed signals as show in Figure 5, where mixed signals are in the left and source signals are in the right.Figure 5After using FastICA, we recover the source, and we show the estimated signals in the left with scale and original signals in the right in Figure6.Figure 6


## 5.1. Example  1

In this example,k=1 and the data set used consists of different realizations of independent signals, with distributions shown in Table 2. Note that this is a large-scale and substantially difficult separation problem, since it involves a Gaussian, various super- and sub-Gaussian symmetric PDFs, as well as asymmetric distributions. In all cases, the number of data samples has also been designed to be relatively small; for example, N=250. The source signals are mixed (noiselessly) with randomly generated full-rank mixing matrices A. The FastICA method is implemented in the so-called simultaneous separation mode whereas the stopping criterion is set to ε=10-4. FastICA is executed using the flexible MGΓD model is used to model the distribution of the unknown sources, while (3.3)–(3.5) are employed to adaptively calculate the necessary parameters of the MGΓD-based score function defined in (4.4).Table 2SourcesMGΓD(α,γ)PDF plotGaussian0.5, 2.0Laplace1.0, 1.0Uniform0.1, 10Gamma0.5, 1.0Now, to show the performance in this case, we consider three source signals (source), where these signals are generated randomly from Weilbull, gamma, and exponential distributions as follows:(5.1)s1=wblrnd(1:250,1),s2=exprnd(1,1,250),s3=gamrnd(1,1,250).
Let the mixing matrix A and unmixing matrix W be defined as follows:(5.2)A=(.56.79.37.75.65.86.17.32.48),W=(0.2128-0.21770.25660.2184-0.19240.16461.9485-1.87981.8687).
By using the equation x=As, we obtain mixed signals as show in Figure 3, where mixing signals are in the left and source signals are in the right.Figure 3After using FastICA, we recover the sources, and we show the estimated signals in the left and original signals in the right in Figure4 with different in scale only.Figure 4


## 5.2. Example  2

In this example,k=2 in which the data set used consists of different realizations of independent signals. Note that each signal not only a Gaussian, super-, and sub-Gaussian PDFs, but it is mixed of this PDFs as shown in Table 3.Table 3SourcesMGΓD(α1,γ1,α2,γ2)PDF plotGamma_Lapalce0.5, 1.0, 1.0, 1.0Gaussian_Laplace0.5, 2.0, 1.0, 1.0Weilbull_Gamma2.0, 4.0, 0.5, 1.0Uniform_Laplace0.1, 10, 1.0, 1.0To show the performance in this case, we consider three source signals (source) where these signals are generated randomly from Gamma_Lapalce (s1), Weilbull (s2), and Gaussian_Laplace (s3)  distributions. Let the mixing matrix A and unmixing matrix W be defined as follows:(5.3)A=(.56.79.37.75.65.86.17.32.48),W=(-0.12932.6579-4.2322-1.4418-3.15339.9930-2.33108.4107-11.6627).
By using the equation X=AS,we obtain mixed signals as show in Figure 5, where mixed signals are in the left and source signals are in the right.Figure 5After using FastICA, we recover the source, and we show the estimated signals in the left with scale and original signals in the right in Figure6.Figure 6


## 6. Algorithm Performance

The separation performance for ICA algorithm is evaluated with the crosstalk error measure(6.1)PI=∑i=1n(∑j=1n|pij|2max⁡1≤l≤n|pil|2-1)+∑j=1n(∑i=1n|pij|2max⁡1≤l≤n⁡|plj|2-1).
Note that here, pij represents the elements of the permutation matrix P=WA, which after assuming that all sources have been successfully separated should ideally reduce to a permuted and scaled version of the identity matrix. The separation performance for the first example is PI=-12.68dB and for second example is PI=-16.68dB.


## 7. Conclusions

We have derived a novel parametric family of flexible score functions, based exclusively on the MGΓD model. To calculate the parameters of these functions in an adaptive BSS setup, we have chosen to maximize the ML equation with the NM optimization method. This alleviates excessive computational cost requirements and allows for a fast practical implementation of the FastICA. Simulation results show that the proposed approach is capable of separating mixtures of signals.



---

*Source: 150294-2010-11-30.xml*