Nongaussian density estimation for the classification of acoustic feature vectors in speech recognition

[0002] 1. Field of the Invention

[0003] The present invention generally relates to speech recognition systems and, more particularly, to the use of EM type algorithms for the estimation of parameters for a mixture model of nongaussian densities. The present invention was motivated by two objectives. The first was to study maximum likelihood density estimation methods for high dimensional data, and the second was the application of the techniques developed in large vocabulary continuous parameter speech recognition.

[0004] 2. Background Description

[0005] Speech recognition systems requires modeling the probability density of feature vectors in the acoustic space of phonetic units. Purely gaussian densities have been know to be inadequate for this purpose due to the heavy tailed distributions observed by speech feature vectors. Se,, for example, Frederick Jelenik, Statistical Methods for Speech Recognition, MIT Press, 1997. As an intended remedy to this problem, practically all speech recognition systems attempt modeling by using a mixture model with gaussian densities for mixture components. Variants of the standard K-mean clustering algorithm are used for this purpose. The classical version (as described by John Hartigan in Clustering Algorithms, John Wiley & Sonse, 1975, and Anil Jain and Richard Dubes in Algorithms for Clustering Data, Prentice Hall, 1988) of the K-means algorithm can also be viewed as an special case of the EM algorithm (as described by A. P. Dempster, N. M. Laird and D. B. Baum in “Maximum likelihood from incomplete data via the EM algorithm”, Journal of Royal Statistical Soc., Ser. B, vol. 39, pp. 1-38, 1997) in the limiting case of gaussian density estimation with variance zero. See, for example, Christopher M. Bishop, Neutral Networks for Pattern Recognition, Cambridge University Press, 1997, and F. Marroquin and J. Girosi, “Some extensions of the K-means algorithm for image segmentation and pattern classification”, MIT Artificial Intelligence Lab. A. I. Memorandum no.1390, January 1993. The only known attempt to model the phonetic units in speech with nongaussian mixture densities is described by H. Ney and A. Noll in “Phoneme modeling using continuous mixture densities”, Proceedings of IEEE Int. Conf. on Acoustics Speech and Signal Processing, pp. 437-440, 1988, where laplacian densities were used in a heuristic based estimation algorithm.

[0006] It is therefore an object of this invention to provide a new statistical modeling paradigm for automatic machine recognition of speech.

[0007] According to this invention, novel mixtures of nongaussian statistical probability densities for modeling speech subword units (e.g., phonemes) and further subcategories thereof, including the transition and output probabilities, are used in a Hidden Markov generative model of speech. We model speech data by building probability densities from functions of the form exp(−t^α/2) for t≧0, α>0. In this notation, the case α=2 corresponds to the gaussian density, whereas the laplacian case considered in Ney et al. corresponds to α=1.

[0008] Furthermore, we focus on a set of four different types of mixture components constructed from a different univariate function. For each of them, a mixture model is then used for a maximum likelihood model of speech data. It turns out that our iterative algorithm can be used for a range of values of α (opposed to fixed α=1 in Ney et al. or α=2 in standard speech recognition systems).

[0009] Our observation is that the distribution of speech feature vectors in the acoustic space are better modeled by mixture models with nongaussian mixture components. In particular, for speech α<1 seems more appropriate, see FIG. 1. To wit, very similar distributions have been noted for the distribution of image gradients by Stuart Geman in “Three lectures on image understanding”, The Center For Imaging Science, Washington State University, video tabe, Sep. 10-12, 1997, and also by David Mumford in pattern theory, Lecture in Directors Series, IBM Yorktown Heights, Feb. 23, 1998.

[0010] The second point to be made is that from a practical standpoint estimation of densities in speech data is accomplished by all the difficulties characteristic of high dimensional density estimation problems. See David W. Scott, Multivariate Density Extimation, Wiley Interscience, 1992, and James R. Thompson and Richard A. Tapia, Nonparametric Function Estimation, modeling and simulation, SIAM Publications, 1997. Feature vectors of dimension fifty or more are typical in speech recognition systems, and consequently, the data can be considered to be highly sparse. In contrast, the literature (see Scott, supra) on multivariate density estimation puts high emphasis on “exploratory data analysis” the goal of which is to glean insight about the densities via visualization of the data. This is not feasible for dimensions of the order of fifty or more, even when projections on lower dimensional spaces are considered.

[0011] The classification/training step for speech recognition which we use can be cast in a general framework. We begin with a parametric family p(x/λ), x ε R^d, λ ε Ω ⊂R^qof probability densities on R^dwith parameters in the manifold Ω in R^q. The method of classification used here begins with k finite subsets of T₁, T₂, . . . , T_kof R^dand consider the problem of deciding which of these sets given vector x ε R^dlies. The method that we employ picks k probability densities p₁=p(.|θ₁), . . . , p_k=p (.|θ_k) from our family and associates with the subset T_lthe probability density p₁, l=1, 2, . . . , k. Then, we say x ε R^dbelongs to the subset T_rif r is the least integer in the set {1, 2, . . . , k} such that $p_r\left(x\right)\ge \max \left\{p_l\left(x\right)\text{:}1\le l\le k\right\}.$

[0012] To use this method we need to solve the problem of determining, for a given finite subset T⊂R^d, a probability density p(.|θ) in our family. This is accomplished by maximum likelihood estimation (MLE). Thus, the likelihood function, for the data T is given by $L\left(\lambda |T\right)=\prod _{y\in T}p\left(y|\lambda \right),\lambda \in \Omega ,$

[0013] and a vector, θ ε Ω, is chosen which maximizes this function over all λ ε Ω (if possible). Generally, a maximum does not exist, (see V. N. Vapnik, The Nature of Statistical Learning Theory, Springer Verlag, 1995, p. 24), and thus, typically as iterative method is used to find a stationary point of the likelihood function. As we shall see, the iteration we use takes the form of a variation of the EM algorithm described by Dempster et al., supra, and Richard A. Redner and Homer Walker, “Mixture densities, maximum likelihood and the EM algorithm”, SIAM Review vol. 26, no. 2, April 1984.

[0014] The foregoing and other objects, aspects and advantages will be better understood from the following detailed description of a preferred embodiment of the invention with reference to the drawings, in which:

[0015] FIG. 1 is a graph showing an optimum mixture component in logarithmic scale;

[0016] FIG. 2 is a graph showing a one-dimensional density profile for different values of α;

[0017] FIGS. 3A, 3B, 3C, and 3D are graphs of level curves of two-dimensional tensor product α-density for different values of α;

[0018] FIG. 4 is a graph of a clipped gaussian curve;

[0019] FIG. 5 is a graph of a gaussian curve with double exponential tails;

[0020] FIG. 6 is a graph demonstrating the four fixed points of the map V_3/2;

[0021] FIG. 7 is a graph of W_3/2(μ) and its minimum as a function of μ;

[0022] FIG. 8 is a graph of W_1/2(μ) and its minimum as a function of μ;

[0023] FIG. 9 is a graph demonstrating the seven fixed points of the map V_½;

[0024] FIG. 10 is a graph showing percent word error for different values of α; and

[0025] FIG. 11 is a flow diagram illustrating a computer implementation of the process according to the invention.

[0026] TAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION We begin with a nonnegative function f defined on the positive real axis R₊ such that all of its moments $\begin{array}{cc}{m}_{\beta }={\int }_{R_{+}}{t}^{\beta }f\left(t\right)dt,& \left(1\right)\end{array}$

[0027] are finite for β≧−½. Here, R₊ denotes the positive real axis. We also use PD^dto denote the space of all positive definite matrices of size d.

[0028] Lemma 1 We define constants ${\rho }_d=\frac{\Gamma \left(\frac{d}{2}\right)}{{\pi }^{\frac{d}{2}}}\frac{{\left(m_{\frac{d}{2}}\right)}^{\frac{d}{2}}}{{\left(m_{\frac{d}{2}-1}\right)}^{\frac{d}{2}+1}}$

[0029] and ${\gamma }_d=\frac{m_{\frac{d}{2}}}{m_{\frac{d}{2}-1}}.$

[0030] Then for any vector μ in R^dand any positive definite matrix Σ of size d, the function $\begin{array}{cc}p\left(x|\mu ,\Sigma \right)={\rho }_d\frac{1}{\sqrt{\det \Sigma }}f\left({{\gamma }_d\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right)\right),x\in R^d& \left(2\right)\end{array}$

[0031] is a probability density function with mean μ and covariance Σ.

[0032] Proof: The proof uses the following identity which holds for any non-negative integer k and positive constant δ ${\int }_{R^d}{\Vert x\Vert }^{k}f\left(\delta {\Vert x\Vert }^2\right)dx=\frac{{\pi }^{\frac{d}{2}}}{\Gamma \left(\frac{d}{2}\right)}{\delta }^{-\frac{k+d}{2}}{m}_{\frac{k+d-2}{2}},$

[0033] where ${\Vert x\Vert }^2=\sum _{l=1}^d{x}_l^2,x=\left(x_1,x_2,\dots ,x_d\right).$

[0034] the verification of this equation uses the fact that volume of the surface of unit sphere in R^dis $\frac{2{\pi }^{\frac{d}{2}}}{\Gamma \left(\frac{d}{2}\right)}.$

[0035] Additionally, we may combine the method of tensor product with this lemma and construct other densities which we shall employ later. The proof of the next lemma is even simpler than the one above.

[0036] Lemma 2 Let x=(x¹, x²), where xⁱε R^d, i=1,2. Suppose p_i(xⁱ|Σⁱ, μⁱ), xⁱε R^d; i=1,2 have the form (2) above. Then $p(\left(x^1,x^2(|\mu ,\Sigma \right)=p_1\left(x^1|{\mu }^1,{\Sigma }^1\right)p_2\left(x^2|{\mu }^2,{\Sigma }^2\right),x^1\in R^{d_1},x^2\in R^{d_2}$

[0037] has mean μ=(μ¹, μ²) and covariance $\Sigma =\left(\begin{array}{cc}{\Sigma }^1& 0\\ 0& {\Sigma }^2\end{array}\right).$

[0038] A special case of this construction is the density $p\left(x|\mu ,\sigma \right)={\rho }_1^d\frac{1}{\sqrt{\prod _{i=1}^d{\sigma }_i}}\prod _{i=1}^{d}f\left({\gamma }_1\frac{{\left(x_i-{\mu }_i\right)}^2}{{\sigma }_i}\right)$

[0039] defined for x=(x₁, x₂, . . . , x_d) ε R^d, where μ=(μ₁, μ₂, . . . , μ_d) and σ=(σ₁, σ₂, . . . , σ_d) are the mean and diagonal covariance, respectively.

[0040] Let us note some special cases.

[0041] For this example, we choose any α ε R₊, define the function $\begin{array}{cc}f\left(t\right)=\exp \left(-t^{\alpha /2}\right),t\in R_{+}& \left(3\right)\end{array}$

[0042] and note for β>−1 that $\begin{array}{cc}{m}_{\beta }=2{\alpha }^{-1}\Gamma \left(\frac{2\beta +2}{\alpha }\right).& \left(4\right)\end{array}$

[0043] Consequently, the function $\begin{array}{cc}{S}_{\beta }\left(x|\mu ,\Sigma \right)=\frac{{\rho }_d\left(\alpha \right)}{\sqrt{\det \Sigma }}\exp \left(-{\left({\gamma }_d\left(\alpha \right)\left(x-\mu \right)\right)}^{\frac{\alpha }{2}}\right),x\in R^d& \left(5\right)\end{array}$

[0044] where ${\rho }_d\left(\alpha \right)=\frac{\alpha }{2}\frac{\Gamma \left(\frac{d}{2}\right)}{{\pi }^{\frac{d}{2}}}\frac{{\Gamma \left(\frac{d+2}{\alpha }\right)}^{\frac{d}{2}}}{{\Gamma \left(\frac{d}{\alpha }\right)}^{\frac{d}{2}+1}}$

[0045] and ${\gamma }_d\left(\alpha \right)=\frac{\Gamma \left(\frac{d+2}{\alpha }\right)}{\Gamma \left(\frac{d}{\alpha }\right)}$

[0046] is a probability density with mean μ and covariance Σ for any α>0. We refer to this as a spherical α-density.

[0047] As a special case of this construction, we choose d=1 and consider the univariate density

[0048] We use this univariate density and construct a multivariate density by the formula

[0049] where x=(x₁, x₂, . . . , x_d) ε R^d. We refer to this as a tensor product α- $\begin{array}{cc}\frac{1}{\sqrt{\sigma }}{\rho }_1\left(\alpha \right)\exp \left(-{\gamma }_1\left(\alpha \right)\left(\frac{{\left(t-\mu \right)}^{\alpha }}{\sigma }\right)\right),t\in R.& \left(6\right)\\ T_{\alpha }\left(x|\mu ,\sigma \right)=\frac{1}{\sqrt{\det \mathrm{diag}\sigma }}{{\rho }_1\left(\alpha \right)}^d\exp \left(-{{\gamma }_1\left(\alpha \right)}^{\frac{\alpha }{2}}\sum _{l=1}^d{\sigma }_l^{-\alpha /2}{\left|x_l-{\mu }_l\right|}^{\alpha }\right),& \end{array}$

[0050] density.

[0051] It is interesting to note the limiting form of these distributions when α→^∞. (We refer to K. Fukunaga, Statistical Pattern Recognition, second edition, Academic Press, 1990, for a related discussion). Since lim_x−0xΓ(x)=1, we conclude that $\underset{\alpha \to \infty }{\lim }{\gamma }_d\left(\alpha \right)=\frac{d}{d+2}$

[0052] and also $\underset{\alpha \to \infty }{\lim }{\rho }_d\left(\alpha \right){\left\{\frac{d}{d+2)\pi }\right\}}^{\frac{d}{2}}\Gamma \left(\frac{d}{2}+1\right).$

[0053] Hence, we obtain for almost all x ε R^dthat $\underset{\alpha \to \infty }{\lim }{S}_{\alpha }\left(x|\mu ,\Sigma \right)=S_{\infty }\left(x|\mu ,\Sigma \right)$

[0054] where $S_{\infty }\left(x|\mu ,\Sigma \right):=\frac{1}{\sqrt{\det \Sigma }}{f}_{\infty }\left({\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right)\right),x\in R^d$

[0055] and for t ε R $f_{\infty }\left(t\right):={\left\{\frac{d}{\left(d+2\right)\pi }\right\}}^{d/2}\Gamma \left(\frac{d+2}{2}\right)\{\begin{array}{cc}1,& \left|t\right|<\frac{d=2}{2}\\ 0,& \mathrm{otherwise}\end{array}$

[0056] Similarly, for the tensor product construction, we have for any x ε R^dthat $\underset{\alpha ->\infty }{\lim }{T}_{\alpha }(x|\mu ,\sigma )=\frac{1}{\sqrt{\prod _{i=1}^d{\sigma }_i}}{g}_{\infty }\left({\Vert \left(x-\mu \right)/\sqrt{\sigma }\Vert }_{\infty }\right)$

[0057] where ${\Vert \left(x-\mu \right)/\sqrt{\sigma }\Vert }_{\infty }=\underset{1\le l\le d}{\max }\left|x_l-{\mu }_l\right|/\sqrt{{\sigma }_l}$

[0058] and for t ε R $g_{\infty }\left(t\right)=\frac{1}{3^{d/2}}\{\begin{array}{cc}1,& \left|t\right|\le \sqrt{3}\\ 0,& \mathrm{otherwise}\end{array}$

[0059] The behavior of S_α and T_α where α→0⁺ are equally interesting. Both approach a “delta” function concentrated at their mean. The precise result we record for the spherical case states that $\begin{array}{cc}\underset{\alpha ->0^{+}}{\lim }{\int }_{R_d}f\left(x\right)S_{\alpha }\left(x\right|\mu ,\Sigma )=f\left(\mu \right)& \left(7\right)\end{array}$

[0060] whenever f is continuous on the R^dwith the property that there exist constants M,N>0 and β ε[0,2], such that for ∥x∥>N $\left|f\left(x\right)\right|\le M{\Vert x\Vert }^{\beta .}$

[0061] Observe that the value β=2 must be excluded since for the function $f\left(x\right)={\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right),x\in R^d$

[0062] the left hand side of equation (7) is one, while the right hand side is zero. The verification of equation (7) follows from a computation of all the moments of S_α. Specifically, we have the following fact.

[0063] Lemma 3 For any β ε (0,2) and α ε R₊, there exist M,N>0 and γ ε (0,1)

[0064] such that ${\int }_{R_d}{\left({\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right)\right)}^{\beta /2}{S}_{\alpha }\left(x\right|\mu ,\Sigma )dx\le M{\lambda }^{1/\alpha }$

[0065] whenever 0<α<N.

[0066] Proof: We proceed as in Lemma 1 and observe for any β>0 that ${\int }_{R_d}{\left({\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right)\right)}^{\beta /2}{S}_{\alpha }\left(x\right|\mu ,\Sigma )dx=\frac{\Gamma \left(\frac{\beta +d}{\alpha }\right)}{{\Gamma \left(\frac{d+2}{\alpha }\right)}^{\beta /2}{\Gamma \left(\frac{d}{\alpha }\right)}^{1-\frac{\beta }{2}}}$

[0067] An application of Stirling's formula (see M. Abramowitz and I. Stegan, Handbook of Mathematical Functions, Dover Publications, New York, Ninth Dover printing, 1972, pp. 254, 6.1.38) $1<\Gamma \left(t+1\right)/\sqrt{2\pi }{t}^{t+½}\exp \left(-t\right)<\exp \left(1/12t\right),t\in R$

[0068] proves the result with $\lambda :=\frac{{\left(\beta +d\right)}^{\beta +d}}{{\left(d+2\right)}^{\beta \left(d+2\right)/2}{d}^{d\left(1-\beta /2\right)}}$

[0069] Note that using the same methods as above, we may conclude that $\underset{\alpha ->0^{+}}{\lim }{\alpha }^2{\gamma }_{\alpha }^d\left(\alpha \right)=\frac{{\left(d+2\right)}^{d+2}}{e{d}^d}$

[0070] and ${\tilde{\rho }}_d\left(\alpha \right)<{\rho }_d\left(\alpha \right)<\exp \left(\alpha /12\right){\tilde{\rho }}_d\left(\alpha \right)$

[0071] where ${\tilde{\rho }}_d\left(\alpha \right):={\alpha }^{-½}\frac{\Gamma \left(\frac{d}{2}\right)}{{\pi }^{\frac{d}{2}+1}}{\left[\frac{d^{1+\frac{d}{2}}}{2{\left(d+2\right)}^{\frac{d}{2}}}\right]}^{\frac{3}{2}}{\left[\frac{\left(d+2\right)}{d}\right]}^{\frac{d\left(d+2\right)}{2\alpha }}$

[0072] The statement for the tensor product case is exactly the same and is proved in the same way. The univariate α-densities with zero mean and unit variance a graphically shown in FIG. 2. Note that the covariance matrix of the tensor product α-density is the diagonal matrix diag σ. Also, when α is equal to 2 and Σ equals diag σ, then both the spherical α-density and the tensor product α-density are identical with the gaussian density. However, it must be noted that the tensor product α-density is different from the spherical α-density for α≠1 even when Σ is diagonal. In fact, assuming zero means and unit covariances in equations (5) and 6), the contours of constant density for the spherical α-density in equation (5) is given

[0073] by $\sum _{i=1}^d{x}_i^2$

[0074] =constant, whereas contours of constant density for the tensor product α-density in equation (6) is given by $\sum _{i=1}^d{x}_i^{\alpha }$

[0075] =constant. The former set of contors are plotted for different values of α in FIGS. 3A to 3D thereby illustrating the difference between the two densities even for diagonal covariances.

[0076] The next example we shall consider is a clipped gaussian density.

[0077] For very ε>0 and t nonegative, we set

[0078] h(t)=max(0,t−ε).

[0079] In this case, for every integer n, the moments of the function f=exp(−h) required for the spherical construction are given by the formula $m_{\beta }=\frac{{\epsilon }^{\beta }}{\beta +1}+\exp \left(\epsilon \right){\int }_{\epsilon }^{\infty }{t}^{\beta }\exp \left(-t\right)dt,\beta >-1,\epsilon >0.$

[0080] Moreover, for $r=\left[\frac{d}{2}\right]\left(i.e.,\mathrm{the}\mathrm{greatest}\mathrm{integer}\le \frac{d}{2}\right)$

[0081] since

[0082] we get that

[0083] where we use the definition (x)₁=Γ(x+1)/Γ(x−1+l). For an example of ${\int }_{\epsilon }^{\infty }{t}^{d/2}\exp \left(-t\right)dt=\{\begin{array}{c}\sum _{l=0}^r{\left(r\right)}_l{\epsilon }^{r-l}\\ \sum _{l=0}^r{\left(\frac{d}{2}\right)}_l{\epsilon }^{\frac{d}{2}-l}+{\left(\frac{d}{2}\right)}_{r+1}\sqrt{\pi }\exp \left(\epsilon \right)\left(1-\mathrm{erf}\left(\sqrt{\epsilon }\right)\right)\end{array}{m}_{\frac{d}{2}}=\{\begin{array}{c}\sum _{l=1}^r{\left(r\right)}_l{\epsilon }^{r-l}\\ \sum _{l=1}^r{\left(\frac{d}{2}\right)}_l{\epsilon }^{\frac{d}{2}-l}+{\left(\frac{d}{2}\right)}_{r+1}\sqrt{\pi }\exp \left(\epsilon \right)\left(1-\mathrm{erf}\left(\sqrt{\epsilon }\right)\right)\end{array}$

[0084] this density see FIG. 4.

[0085] The final example we shall consider is a gaussian density with a double exponential tails introduced by P. Huber (see for instance D. M. Titterington, A. F. M. Smith and U. E. Markov, Statistical Analysis of Finite Mixture Distributions, Wiley Interscience, New York, 1985, p.23).

[0086] For every ε>0 and t nonegative, we set $h\left(t\right)=\{\begin{array}{cc}t,& t\le \epsilon \\ 2\sqrt{\epsilon t}-\epsilon ,& t\le \epsilon \end{array}$

[0087] The required moments are given by the formula $m_{\frac{d}{2}}=\Gamma \left(\frac{d}{2}+1\right)-{\int }_{\epsilon }^{\infty }{t}^{\frac{d}{2}}\exp \left(-t\right)dt+\frac{1}{2^{d+1}{\epsilon }^{\frac{d}{2}}}\sum _{l=0}^d{\left(d\right)}_l{\epsilon }^{d-l}.$

[0088] The corresponding density is a gaussian whose tails have been adjusted to have an exponential decay, see FIG. 5 for its plot.

[0089] Maximum Likelihood Estimation

[0090] Let T={x¹, x², . . . , x^N} R^dbe independent samples from a probability density function p. Then the likelihood function of the data is $\prod _{k=1}^{N}p\left(x^k\right).$

[0091] Given a function f as above, we discuss the problem of determining the value of the mean μ and R^dand a covariance matrix Σ which maximizes the likelihood function for the density $p\left(x\right|\mu ,\Sigma )={\rho }_d\frac{1}{\sqrt{\det \Sigma }}f\left({{\gamma }_d\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right)\right),t\in R^d.$

[0092] It is convenient to express f in the form

[0093] f(t)=exp(−h(t)), t ε R

[0094] and consider the negative of the log-likelihood function $\begin{array}{cc}-\log L=-N\log {\rho }_d-\frac{N}{2}\log \det {\Sigma }^{-1}+\sum _{k=1}^N\left({{\gamma }_d\left(x^k-\mu \right)}^t{\Sigma }^{-1}\left(x^k-\mu \right)\right)& \left(8\right)\end{array}$

[0095] where x¹, x², . . . , x^Nare the given data. We introduce the positive definite matrix Γ, $\Gamma ={\gamma }_d^{½}{\Sigma }^{-½},$

[0096] the vector

[0097] and function F

[0098] Therefore, we can express the log-likelihood in the alternate form $\varphi ={\gamma }_d^{½}\Gamma \mu$ $F:=F\left(\Gamma ,\varphi \right):=-N\log \det \Gamma +\sum _{k=1}^{N}h\left({\Vert {\Gamma x}^k-\varphi \Vert }^2\right).-\log L=-N\log {\rho }_d+\frac{N}{2}\log {\gamma }_d+F.$

[0099] Thus, maximizing L over μ ε R^dand Σ ε PD^dis equivalent to minimizing F over Γ ε PD^dand φ ε R^d. Our first lemma identifies a class of univariate functions f, for which L achieves its maximum R^d×PD^d. To this end, it is convenient to introduce the linear subspace of R^ddefined by $\begin{array}{cc}M\left(T\right):=\left\{\sum _{k=1}^N{c}_k{x}^k:\sum _{k=1}^N{c}_k=0\right\}.& \left(9\right)\end{array}$

[0100] Lemma 4 Suppose there exist constants a, b and c, with b, c positive such that $h\left(t\right)\ge a+b{t}^c,t\in R_{+}$

[0101] and

[0102] dim M=d

[0103] Then L has a maximum on R^d×PD^d.

[0104] Proof: Since F is a continuous function of Γ and Φ (and F is ∞ when Γ is singular), it suffices to prove that F goes to infinity as (Γ,Φ)→∞. To this end we observe that the function

[0105] vanishes if and only if both Γ and μ are zero. We conclude that there exists $\sum _{k=1}^N{\Vert {\Gamma x}^k-\varphi \Vert }^c$

[0106] a positive constant b′ such that $F\ge -N\log \det \Gamma +aN+b^{\prime }\left({\Vert \Gamma \Vert }^c+{\Vert \varphi \Vert }^c\right)$

[0107] where ∥Γ∥ denotes the Frobenius norm of Γ. By the Hadamard inequality (see R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press,1985, p. 477), we obtain that det Γ≦∥Γ∥^d, and so $F\ge -dN\log \Vert \Gamma \Vert +aN+b^{\prime }\left({\Vert \Gamma \Vert }^c+{\Vert \varphi \Vert }^c\right)$

[0108] which proves the result.

[0109] The next lemma identifies a condition on h so that L has a unique maximum.

[0110] Lemma 5 Let H(t) :=h(t²), t ε R₊ and suppose that H is convex on R₊. Then L has at most one maximum.

[0111] Proof: Recall that as a function of Γ ε PD^d, det Γ is strictly log concave (see R. A. Horn et al., supra). Hence, the function F(Γ,φ) defined above is a strictly convex function of Γ and φ. This proves the result.

[0112] When h(t)=t, the corresponding spherical density is a gaussian distribution, and in this case, the MLE of Σ and μ are given by the formula $\mu =\frac{1}{N}\sum _{k=1}^N{x}^k;\Sigma =\frac{1}{N}\sum _{k=1}^N\left(x^k-\mu \right){\left(x^k-\mu \right)}^t.$

[0113] Except for this case, the optimal value of Σ and μ must be found by some iterative method.

[0114] We begin our discussion of the problem with the univariate case of α-density, a case of practical importance to us here. Thus, we have scalar data {x₁, x₂, . . . , x_N}, and the negative log-likelihood is $-\log L=N\log {\rho }_1\left(\alpha \right)+\frac{N}{2}\log \sigma +{{\gamma }_1\left(\alpha \right)}^{\alpha /2}{\sigma }^{-\alpha /2}\sum _{k=1}^N{\left|x_k-\mu \right|}^{\alpha }.$

[0115] Notice here that the minimum over μ ε R and σ ε R₊ can be performed separately. For this reason, we set $\omega ={\min }_{\mu }{\left(\frac{1}{N}\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha }\right)}^{\frac{1}{\alpha }}={\left(\frac{1}{N}\sum _{i=1}^N{\left|x_i-\hat{\mu }\right|}^{\alpha }\right)}^{\frac{1}{\alpha }}.$

[0116] Since the derivative of (−log L) with respect to σ is given by (at μ={circumflex over (μ)}) $\frac{d}{d\sigma }\left(-\log L\right)=\frac{N}{2\sigma }\left(1-\frac{\alpha {{\gamma }_1\left(\alpha \right)}^{\frac{\alpha }{2}}{\omega }^{\alpha }}{{\sigma }^{\frac{\alpha }{2}}}\right)$

[0117] the unique value of Σ at which the minimum occurs is $\begin{array}{cc}\hat{\sigma }={\alpha }^{\frac{2}{\alpha }}{\gamma }_1\left(\alpha \right){\omega }^2.& \left(10\right)\end{array}$

[0118] Note that the hypothesis of Lemma 5 is satisfied for α≧1. Since in this case {circumflex over (μ)} is unique, our computation confirms directly that −log L has a unique minimum when α≧1 (but not for 0<α<1, see below). Besides the case α=2, corresponding to the gaussian distribution, we note that when α=1 we have ρ₁(α)=γ₁(α)=2 and so $\sigma =2{\left(\frac{1}{N}\sum _{i=1}^N\left|x_i-\tilde{\mu }\right|\right)}^2$

[0119] where {circumflex over (μ)} is the median of the data.

[0120] In the general case, the minimum value of −log L has the form $-\log L=N\left\{\kappa +\log \omega \right\}$

[0121] where $\kappa :={\alpha }^{-1}+\left(1-{\alpha }^{-1}\right)\log \alpha +\log \Gamma \left({\alpha }^{-1}\right)$

[0122] is a constant independent of the data. Therefore, a comparison of the likelihoods of two data sets is equivalent to a comparison of their corresponding w values.

[0123] There remains the problem of determining {circumflex over (μ)}. Because of its role in the EM algorithm, see below, we accomplish this by the update formula $\begin{array}{cc}\hat{\mu }=\frac{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha -2}{x}_i}{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha -2}}.& \left(11\right)\end{array}$

[0124] For α<2, care must be taken in the use of this formula. If μ=x_lfor some l, l=1, . . . , N, then {circumflex over (μ)}=x_l. Otherwise, division by zero does not occur in equation (11). The method always converges to the unique minimum ũ of W independent of the choice of the initial value used when α≧2. the reason for this is two-fold: first, |{circumflex over (μ)}|≦max{|x_i|; i=1,2, . . . , N} and secondly, the function $\begin{array}{cc}{W}_{\alpha }\left(\mu \right):=\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha },\mu \in R& \left(12\right)\end{array}$

[0125] is a strictly convex function of μ. Thus, {circumflex over (μ)}=μ in equation (11) if and only if the derivative of the function of μ in equation (12) vanishes at {circumflex over (μ)}. For α<2, the issue of convergence of the iteration is delicate. We divide our comments into two cases. The first case we consider is 1<α2. The iteration will converge for a generically chosen initial data. We illustrate this phenomenon with an example. When α=3/2, N=3, x₁=0, x₂=¼, x₃=1, the function in equation (12) shown in FIG. 7 has its unique minimum at =0.32. As expected, this value of μ is a fixed point of the iteration of equation (11). In fact, FIGS. 6 and 7 demonstrate that there is a value of μ close t 0.1 so that the updated value {circumflex over (μ)} is equal to ¼. Since this is also a fixed point of the mapping in equation (11), all further updates does not change the value of μ. Moreover, starting at any value of μ except for the (countable) pre-images of x₂(which has zero as its only accumulation points) the iteration of equation (11) will converge geometrically fast to the minimum of equation (12).

[0126] We do not present a complete analysis of the convergence of this method here. The following observation provides some insight into what might be expected.

[0127] Lemma 6 Suppose 1<α<2, and let {μ_n: n=1, 2, . . . } be a sequence generated by the iteration of equation (11), which converges to some μ. Either μ is equal the unique minimum of the function appearing in equation (12) or there exists an l in the set {1, 2, . . . , N} such that all the elements of the sequence except for a finite number equal x_l.

[0128] Proof: If μis not in the set {x_j: j=1, 2, . . . , N} then μ must be the unique minimum of the function in equation (12). In the case that μ=x₁for some l=1,2, . . . , N, there exists a positive number ε such that whenever |x−x_l|<ε it follows that |V_α(x)−x_l|≦|x−x_l|, where V_α(x) is defined to be $\begin{array}{cc}{V}_{\alpha }\left(x\right)=\frac{\sum _{i=1}^N{\left|x_i-x\right|}^{\alpha -2}{x}_i}{\sum _{i=1}^N{\left|x_i-x\right|}^{\alpha -2}}.& \left(13\right)\end{array}$

[0129] In other words, all the points x₁, x₂, . . . , x_Nare repelling fixed points of the map of equation (13). Assuming that this be the case, it follows that all but a finite number of this sequence must equal x_l. Thus, the proof is completed by demonstrating the existence of the ε mentioned above. To this end, we merely observe that the right and left derivatives of Vα(x) at any data point x_l, l=1,2, . . . , N is infinite.

[0130] Let is not comment on the case 0<α<1. In general, the iteration does converge to the minimum of the function W for most initial values. In this regard, we point out that W is strictly concave between each data point and has infinite left and right derivatives at the data points. Hence, there is a unique local maximum between each data point, but in general, there is more than one minimum of W. This means that, in general, the maximum of the likelihood occurs for more than one value of mean and variance.

[0131] Our computational experience seems to indicate that the iteration of equation (11) will converge to the data point within the between two consecutive local maxima of the function W, in which the initial value lies. We illustrate this with an example when α=0.5 and x_i=i/10, i=1, . . . 10. The function W is shown in FIG. 8 and the function in equation (13) for α=½ is shown in FIG. 9.

[0132] Other iterations besides equation (11) can be used to find {circumflex over (μ)} as well. For example, since the function appearing in equation (12) is convex, one can use a search method combined with Newton iteration to provide a fast globally convergent algorithm, locating the minimum of this function.

[0133] Returning to the general spherical model, we address the question of iterative methods for finding the MLE. To this end, we compute the variation of log L relative to μ ${\nabla }_{\mu }\log L=2{\gamma }_d{\Sigma }^{-1}\left(\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right)\left(x^k-\mu \right)\right)$

[0134] where $Q\left(x\right|\mu ,\Sigma ):={{\gamma }_d\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right),x\in R^d$

[0135] and for the variation relative to Σ we have that ${\nabla }_{{\Sigma }^{-1}}\log L=\sum _{k=1}^N\left\{\frac{1}{2}\Sigma -{\gamma }_d{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right)\right)\left(x^k-\mu \right){\left(x^k-\mu \right)}^t\}.$

[0136] Thus, (μ,Σ) is a stationary point of L is and only if $\begin{array}{cc}\mu =\frac{\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))x^k}{\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))}& \left(14\right)\end{array}$

[0137] and $\begin{array}{cc}\Sigma =\frac{2{\gamma }_d}{N}\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))\left(x^k-\mu \right){\left(x^k-\mu \right)}^t.& \left(15\right)\end{array}$

[0138] These formulas suggest the iteration scheme in which the updates {circumflex over (μ)}and {circumflex over (Σ)} for the mean μ and the variance σ are respectively given by $\hat{\mu }=\frac{\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))x^k}{\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))}\hat{\Sigma }=\frac{2{\gamma }_d}{N}\sum _{k=1}^N{h}^{\prime }\left(Q\left(x^k\right|\mu ,\Sigma \right))\left(x^k-\mu \right){\left(x^k-\mu \right)}^t.$

[0139] which are related to the EM algorithm (see A. P. Dempster, N. M. Laird and d. B. Baum, “Maximum likelihood from incomplete data via the EM algorithm”, Journal of Royal Statistical Soc.,Ser. B, vol. 39, pp.1-38, 1977, and Richard Redner and Homer Walker, “Mixture densities, maximum likelihood and the EM algorithm, SIAM Review, vol. 26, no. 2, April 1984), discussed below.

[0140] Before getting to this matter, we comment on the convergence of this iteration. All our remarks concern the univariate case. In this case, we rewrite the iteration in terms of the parameters μ and σ corresponding to scalar data {x₁, x₂, . . . , x_N}; thus $\begin{array}{cc}\hat{\mu }=\frac{\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right)x^k}{\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right)}& \left(16\right)\end{array}$

[0141] and $\begin{array}{cc}\hat{\sigma }=\frac{2{\gamma }_1}{N}\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right){\left(x_k-\mu \right)}^2& \left(17\right)\end{array}$

[0142] In general, this iteration does not converge. To see this, we specialize it to the case of Example 1 and observe that in this case, it has the form $\begin{array}{cc}\hat{\sigma }=A\left(\mu \right)/{\sigma }^{\beta }& \left(18\right)\end{array}$

[0143] where $A\left(\mu \right):=\frac{{{\mathrm{\alpha \gamma }}_1\left(\alpha \right)}^{\alpha /2}}{N}{W}_{\alpha }\left(\mu \right),\beta :=\alpha /2-1,$

[0144] and $\begin{array}{cc}\hat{\mu }=V_{\alpha }\left(\mu \right).& \left(19\right)\end{array}$

[0145] We have already pointed out that the iteration for the mean μ converges independent of the initial guess (see discussion after equation (11)) when α≧2. To simplify our exposition, we consider the case that the initial value in the iteration for the mean is chosen to be the minimum of the function appearing in equation (12). In this case, the updated value of the mean and hence also the constant appearing in equation (18) does not change. Let us call this positive constant A. Also, for simplicity, we choose α=4. Thus, the iteration equation (18) takes the form It follows directly, that the solution of this iteration is given by the formula ${\sigma }_{n+1}=\frac{A}{{\sigma }_n^2};n=1,2,\dots$ ${\sigma }_{n+1}={A^{\frac{1}{3}}\left(\frac{{\sigma }_1}{A^{\frac{1}{3}}}\right)}^{{\left(-1\right)}^n{2}^n},n=1,2,\dots$

[0146] The value of σ we seek is A^⅓, (see equation (10)). However, this iteration fails to produce this value except when σ₁is chosen to be A^⅓. Of course, in this simple example, it is not recommended to use the update formula (18), since optimum value of a can be directly obtained (see equation (10)). This example does suggest a means to modify the iteration of equation (17) under certain conditions on the function h. We present this observation in the following lemma.

[0147] Suppose the following four conditions hold.

[0148] The hypothesis of Lemma 5 holds.

[0149] There are nonnegative constants a, b, c such that 0<h′(t)≦a+bf; tεR⁺.

[0150] β:=lim_f−0+fh′(t−1)>0

[0151] The data set {x_k: k=1, 2, . . . , N} consists of at least two points. Then there exists a positive constant κ such that, if σ is at most κ, then {circumflex over (σ)} defined by the formulas $\hat{\mu }=\frac{\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right)x_k}{\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right)}$ ${\hat{\sigma }}^{1+c}=\frac{2{\gamma }_1}{N}\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right){\left(x_k-\mu \right)}^2{\sigma }^c$

[0152] is also at most κ. Moreover, this iteration converges whenever the initial value of σ is at most κ.

[0153] Proof: Set M=max{|xi|: i=1, 2, . . . , N} and first observe that |{circumflex over (μ)}|≦M. Similarly, we obtain that ${\hat{\sigma }}^{1+c}\le u+v{\sigma }^c$

[0154] where μ and ν are positive constants depending on a, b, c, M but not on σ or N. Set $\kappa =\max \left\{2v,{\left(\frac{u}{v}\right)}^{\frac{1}{c}}\right\}.$

[0155] We claim that if σ≦κ, then {circumflex over (σ)}≦κ. If {circumflex over (σ)}≦σ, there is nothing to prove. When {circumflex over (σ)}≧σ, then the above inequality implies that ${\hat{\sigma }}^{1+c}\le u+v{\hat{\sigma }}^c.$

[0156] When $\hat{\sigma }\le {\left(\frac{u}{v}\right)}^{\frac{1}{c}},$

[0157] then again there is nothing to prove. For $\hat{\sigma }\ge {\left(\frac{u}{v}\right)}^{\frac{1}{c}},$

[0158] we have that b{circumflex over (σ)}^c≧μ and so $2v\hat{\sigma }\ge u+v{\hat{\sigma }}^c\ge u+v{\sigma }^c\ge {\hat{\sigma }}^{1+c};$

[0159] that is, {circumflex over (σ)}≦2ν. Hence, under all circumstances we conclude that {circumflex over (σ)}≦κ.

[0160] Now, let {(μ_n, σ_n): n=1, 2, . . . , } be a sequence produced by this iteration, where σ₁is at most κ. Then by our remarks so far, both the sequence of means and variances have a subsequence which converge to some limit, σ_opt, respectively. From the variances update formula and conditions three and four, we have that $\begin{array}{c}{\sigma }_{\mathrm{opt}}^{1+v}=\underset{\_}{{\lim }_{n->\infty }}{\sigma }_n^{1+c}\\ =\frac{2{\gamma }_1}{N}{\lim }_{n->\infty }\sum _{k=1}^N{h}^{\prime }\left(\frac{{{\gamma }_1\left(x_k-\mu \right)}^2}{\sigma }\right){\left(x_k-\mu \right)}^2{\sigma }_n\\ =\frac{2{\gamma }_1}{N}\sum _{k=1}^N{\left(x_k-{\mu }_{\mathrm{opt}}\right)}^{2c+2}>0.\end{array}$

[0161] Hence, we conclude that σ_opt>0. This means that μ_opt, σ_optare stationary points of the likelihood function. However, the first condition guarantees the uniqueness of the stationary point. Therefore, the iterate converges as long as σ₁is at most κ.

[0162] In this section we describe the probability densities which we shall use in our classification experiments. They are built up from mixture models based on the four densities described earlier. Constructing mixture models from a given parametric family is a general procedure. Specifically, starting with a parametric family p(·|λ), λεΩ⊂R^qof probability densities on R^d, a mixture model with m mixtures with mixture weights w_i; i=1, 2, . . . , m has the form $P(x|\omega ,\Lambda ):=\sum _{i=1}^m{\omega }_{ip}(x|{\lambda }_i),\sum _{i=1}^m{\omega }_i=1,{\omega }_i\ge 0,$

[0163] where Λ=(λ₁, λ₂, . . . , λ_m) is chosen from Ω^m=Ω×Ω×. . . ×Ω, m times and ω=(ω₁, ω₂, . . . , ω_m) is a vector with nonnegative components which sum to one. The densities p(·|λ₁), p(·|λ₂), . . . , p(·|λ_m) are called the mixture components P. Note the linear dependency of the mixture model P on the mixture weights.

[0164] Given data vectors {x^k: k=1, 2, . . . , N} in R^d, the goal here is to choose the parameter ΛεΩ^mand the mixture weights ω=(ω₁, ω₂, . . . , ω_m) to maximize the log-likelihood function $\begin{array}{cc}\log L=\sum _{k=1}^{N}P(x^k|\omega ,\Lambda ).& \left(20\right)\end{array}$

[0165] To this end, we first compute the variation of the log-likelihood relative to the mixture weight ω_l, l=1, 2, . . . , m which is given by $\frac{\partial \log L}{\partial {\omega }_l}=\sum _{k=1}^N\frac{p(x^k|{\lambda }_l)}{P(x^k|\omega ,\Lambda )},l=1,2,\dots ,m.$

[0166] Keeping in mind that the mixture weights sum to unity, we are led to the equation $\begin{array}{cc}\frac{1}{N}\sum _{k=1}^N\frac{p(x^k|{\lambda }_l)}{P(x^k|\omega ,\Lambda )},l=1,2,\dots ,m.& \left(21\right)\end{array}$

[0167] This equation for the mixture weights holds independent of the form of the mixture components and suggests the importance of the posterior probabilities $\begin{array}{cc}{P}_l(x|\omega ,\Lambda )=\frac{{\omega }_{l}p(x|{\lambda }_l)}{P(x|\omega ,\Lambda )},l=1,2,\dots ,m,x\in R^d& \left(22\right)\end{array}$

[0168] Thus, equation (21) becomes $\begin{array}{cc}\frac{1}{N}\sum _{k=1}^N{P}_l(x|\omega ,\Lambda )={\omega }_l,l=1,2,\dots ,m& \left(23\right)\end{array}$

[0169] and by the definition of the posterior probability we have that $\begin{array}{cc}\sum _{l=1}^m{P}_l(x|\omega ,\Lambda )=1,x\in R^d& \left(24\right)\end{array}$

[0170] To proceed further, we specify the mixture components to be of the spherical type as in equation (2). In this case, $\log p(x|\mu ,\Sigma )=\frac{1}{2}\log \det {\Sigma }^{-1}-h\left(Q(x|\mu ,\Sigma \right))+\log {\rho }_d,x\in R^d$

[0171] where $Q(x|\mu ,\Sigma ):={{\gamma }_d\left(x-\mu \right)}^t{\Sigma }^{-1}\left(x-\mu \right),x\in R^d$

[0172] Hence, we conclude that ${\nabla }_{\mu }\log p(x|\mu ,\Sigma )=2{\gamma }_d{h}^{\prime }(Q(x|\mu ,\Sigma )){\Sigma }^{-1}\left(x-\mu \right),x\in R^d$

[0173] and ${\nabla }_{{\Sigma }^{-1}}\log p(x|\mu ,\Sigma )=\frac{1}{2}\Sigma -{\gamma }_d{h}^{\prime }(Q(x|\mu ,\Sigma ))\left(x-\mu \right){\left(x-\mu \right)}^t,x\in R^d.$

[0174] Thus, for l=1, 2, . . . , m, the stationary equation for the m mixture means are $\begin{array}{cc}{\mu }^l=\frac{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))x^k}{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))},& \left(25\right)\end{array}$

[0175] and for the variance $\begin{array}{cc}{\Sigma }^l=2{\gamma }_d\frac{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))\left(x^k-{\mu }^l\right){\left(x^k-{\mu }^l\right)}^t}{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )}& \left(26\right)\end{array}$

[0176] where Ψ=(μ¹, Σ¹, . . . , μ^m, Σ^m) is the parameter vector of all the means and matrices. These equations suggest the following update formulas for the means $\begin{array}{cc}{\hat{\mu }}^l=\frac{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))x^k}{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))},& \left(27\right)\end{array}$

[0177] and for the variance $\begin{array}{cc}{\hat{\Sigma }}^l=2{\gamma }_d\frac{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )h^{\prime }\left(Q(x^k|{\mu }^l,{\Sigma }^l\right))\left(x^k-{\mu }^l\right){\left(x^k-{\mu }^l\right)}^t}{\sum _{k=1}^N{P}_l(x^k|\omega ,\psi )}& \left(28\right)\end{array}$

[0178] with l=1, 2, . . . , m.

[0179] Since the problem under consideration is an example of an incomplete data estimation problem (see, for example, the recent book of Christopher M. Bishop, Neural Networks for Pattern Recognition,Cambridge University Press, 1997, for details), the general philosophy of the EM algorithm as described by A. P. Dempster et al., supra, applies. We briefly review the details of the EM algorithm here as it applies to our context. In this regard, we find the point of view of L. E. Baum, Ted Petrie, George Soules, and Norman Weiss in their article “A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains”, The annals of Math. Stat.,vol. 41, no. 1, pp. 164-171, 1970, particularly illuminating. We adopt the notation from that article and let X be a totally finite measure space with measure μ. Let Θ be a subset of some euclidian space and for every θεΘ, we suppose that p(·, θ) is a positive real valued function on X. Let $\begin{array}{cc}P\left(\theta \right)={\int }_{X}p\left(x,\theta \right)d\mu \left(x\right)& \left(29\right)\end{array}$

[0180] and $\begin{array}{cc}Q\left(\theta ,\stackrel{\_}{\theta }\right)={\int }_{X}p\left(x,\theta \right)\log p\left(x,\stackrel{\_}{\theta }\right)d\mu \left(x\right).& \left(30\right)\end{array}$

[0181] Using the concavity of the logarithm, it was shown in L. E. Baum, et al., supra, that P(θ)≦P({overscore (θ)}) whenever Q(θ,θ)≦Q(θ,{overscore (θ)}). Let us suppose that we are given N functions of the form $P_l\left(\theta \right)={\int }_X{p}_l\left(x,\theta \right)d\mu \left(x\right),l=1,2,\dots ,N.$

[0182] We form from these functions the generalized likelihood function $L\left(\theta \right)=\prod _{i=1}^N{P}_l\left(\theta \right)$

[0183] and observe that $L\left(\theta \right)={\int }_{X^N}{p}^N\left(x,\theta \right)d{\mu }_N\left(x\right)$