Impulsivity estimates of mixtures of the power exponential distrubutions in speech modeling

[0002] The present invention generally relates to the technology of speech recognition and, more particularly, to a parametric family of multivariate density functions formed by mixture models from univariate functions for modeling acoustic feature vectors used in automatic recognition of speech.

[0003] Most pattern recognition problems require the modeling probability density of feature vectors in feature space. Specifically, in the problem of speech recognition, it is necessary to model the probability density of acoustic feature vectors in the space of phonetic units. Purely Gaussian densities have been known to be inadequate for this purpose due to the heavy tailed distributions observed by speech feature vectors. See, 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-means clustering algorithm are used for this purpose. The classical version the K-means algorithm as described by John Hartigan in Clustering Algorithms, John Wiley & Sons (1975), and Anil Jain and Richard Dubes in Algorithms for Clustering Data, Prentice Hall (1988), can also be viewed as a special case of the expectation-maximization (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, 1997) for mixtures of Gaussians with variances tending to zero. See also Christopher M. Bishop, Neural 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 attempt to model the phonetic units in speech with non-Gaussian 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, 1998, where Laplacian densities were used in a heuristic base estimation algorithm.

[0004] S. Basu and C. A. Micchelli in “Parametric density estimation for the classification of acoustic feature vectors in speech recognition”, Nonlinear Modeling: Advanced Black-Box Techniques(Eds. J. A. K. Suykens and J. Vandewalle), pp. 87-118, Kluwer Academic Publishers, Boston (1998), attempted to model speech data by building probability densities from a given univariate function h(t) for t≧0. Specifically, Basu and Micchelli considered mixture models from component densities of the form $\begin{array}{cc}p\left(x❘u,\sum \right)={\rho }_d\frac{1}{\sqrt{\det \sum }}\exp (-\left(h\left(Q\left(x\right)\right)\right),x\in R^d\mathrm{where}& \left(1\right)\\ Q\left(x\right)={{\gamma }_d\left(x-\mu \right)}^t{\sum }^{-1}\left(x-u\right),x\in R^d,& \left(2\right)\\ m_{\beta }={\int }_{R_{+}}{t}^{\beta }f\left(t\right)dt,& \left(3\right)\end{array}$

[0005] (when the integral is finite and R₊ denotes the positive real axis) $\begin{array}{cc}{\rho }_d=\frac{\Gamma \left(\frac{d}{2}\right){\left(m_{\frac{d}{2}}\right)}^{\frac{d}{2}}}{{{\pi }^{\frac{d}{2}}\left(m_{\frac{d}{2}-1}\right)}^{\frac{d}{2}+1}},\mathrm{and}& \left(4\right)\\ {\gamma }_d=\frac{m_{\frac{d}{2}}}{d{m}_{\frac{d}{2}-1}}.& \left(5\right)\end{array}$

[0006] If the constraints ρ_dand γ_dare positive and finite, then the vector μ∈R^dand the positive definite symmetric d×d matrix Σ are the mean and the covariance of this density. Particular attention was given to the choice h(t)=t^α/2, t>0, α>0; the case α=2 corresponds to the Gaussian density, whereas the Laplacian case considered by H. Ney and A. Noll, supra, corresponds to α=1. Smaller values of α correspond to more peaked distributions (α→0 yields the δ function), whereas larger values of α correspond to distributions with flat tops (α→∞ yields the uniform distribution over elliptical regions). For more details about these issues see S. Basu and C. Micchelli, supra. This particular choice of densities has been studied in the literature and referred to in various ways; e.g., α-stable densities as well as power exponential distributions. See, for example, E. Gòmez, M. A. Gòmez-Villegas, and J. M. Marin, “A multivariate generalization of the power exponential family of distributions”, Comm. Stat.—Theory Meth. 17(3), pp. 589-600, 1998, and Owen Kenny, Douglas Nelson, John Bodenschatz and Heather A. McMonagle, “Separation of nonspontaneous and spontaneous speech”, Proc. ICASSP, 1998.

[0007] In S. Basu and C. Micchelli, supra, an iterative algorithm having the expectation-maximization (EM) flavor for estimating the parameters was obtained and used for a range of fixed values of α (as opposed to the choice of α=1 in H. Ney and A. Noll, supra, and α=2 in standard speech recognition systems). A preliminary conclusion from the study in S. Basu and C. Micchelli was that the distribution of speech feature vectors in the acoustic space are better modeled by mixture models with non-Gaussian mixture components corresponding to α<1. As a consequence of these encouraging results, we became interested in automatically finding the “best” value of α directly from the data. It is this issue that is the subject of the present invention.

[0008] It is therefore an object of the present invention to provide a parametric family of multivariate density functions formed by mixture models from univariate functions of the type exp(−|x|^β) for modeling acoustic feature vectors used in automatic recognition of speech.

[0009] According to the invention, the parameter β is used to measure the non-Gaussian nature of the data. In the practice of the invention, β is estimated from the data using a maximum likelihood criterion. Among other things, there is a balance between β and the number of data points N that must be satisfied for efficient estimation. The computer implemented method for automatic machine recognition of speech iteratively refines parameter estimates of densities comprising mixtures of power exponential distributions whose parameters are means (μ), variances (σ), impulsivity numbers (α) and weights (w). The iterative refining process begins by predetermining initial values of the parameters μ, σ and w. Then, {circumflex over (μ)}^l, {circumflex over (σ)}^lderived from the following equations ${\mu }_i^l=\frac{\sum _{k=1}^N{\left(\sum _{j=1}^d\frac{{\left(x_j^k-{\hat{\mu }}_j^l\right)}^2}{{\hat{\sigma }}_j^l}\right)}^{{\hat{\alpha }}^l/2-1}{A}_{\mathrm{lk}}{x}_i^k}{\sum _{k=1}^N{\left(\sum _{j=1}^d\frac{{\left(x_j^k-{\hat{\mu }}_j^l\right)}^2}{{\hat{\sigma }}_j^l}\right)}^{{\hat{\alpha }}^l/2-1}{A}_{\mathrm{lk}}}\mathrm{and}$ ${\sigma }_i^l=\frac{{\hat{\alpha }}^l{{\gamma }_d\left({\hat{\alpha }}^l\right)}^{{\hat{\alpha }}^l/2}\sum _{k=1}^N{\left(\sum _{j=1}^d\frac{{\left(x_j^k-{\hat{\mu }}_j^l\right)}^2}{{\hat{\sigma }}_j^l}\right)}^{{\hat{\alpha }}^l/2-1}{A_{\mathrm{lk}}\left(x_i^k-{\hat{\mu }}_i^l\right)}^2}{A_l}$

[0010] for i=1, . . . ,d and l=1, . . . ,m. Then σ is updated by assuming that θ=(μ,σ,α), {circumflex over (θ)}=({circumflex over (μ)},{circumflex over (σ)},{circumflex over (α)}) and letting H(μ,σ)=E_{{circumflex over (θ)}}(log f(·|θ)), in which case H has a unique global maximum at μ={circumflex over (μ)}, σ={circumflex over (σ)} where $\beta \left(\alpha ,\hat{\alpha }\right)={\left\{\frac{\alpha \Gamma \left(\frac{\alpha +1}{\hat{\alpha }}\right)}{\Gamma \left(\frac{1}{\hat{\alpha }}\right)}\right\}}^{\frac{2}{\alpha }}\frac{\Gamma \left(\frac{3}{\alpha }\right)\Gamma \left(\frac{1}{\hat{\alpha }}\right)}{\Gamma \left(\frac{3}{\hat{\alpha }}\right)\Gamma \left(\frac{1}{\alpha }\right)}$

[0011] The l dimension is set by μ^l={circumflex over (μ)}^l, σ^l={circumflex over (σ)}^l, and α^l={circumflex over (α)}^l. Finally, the convergence of a log likelihood function B(α) of the parameters is determined in order to get final values of μ, σ and α. The B(α) is $B\left(\Lambda ,\hat{w},\hat{\Lambda }\right)=\sum _{l=1}^m{B}^l\left(\Lambda ,\hat{w},\hat{\Lambda }\right)$

[0012] where $B^l\left(\Lambda ,\hat{w},\hat{\Lambda }\right)=\sum _{k=1}^N{A}_{\mathrm{lk}}\left(\frac{1}{2}\left(\sum _{i=1}^d\log {\sigma }_i^l\right)+\log {\rho }_d\left({\alpha }^l\right)-{\left({\gamma }_d\left({\alpha }^l\right)\right)}^{{\alpha }^l/2}{\left(\sum _{i=1}^d\frac{{\left(x_i^k-{\mu }_i^l\right)}^2}{{\sigma }_i^l}\right)}^{{\alpha }^l/2}\right).$

[0013] The file of this patent contains at least one drawing executed in color. Copies of this patent with color drawing(s) will be provided by the Patent and Trademark Office upon request and payment of the necessary fee.

[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 the plot of average log-likelihood for −1<μ<1 and 0<α<5;

[0016] FIG. 2 is a graph showing the plot of average log-likelihood for −1<μ<1 and 0<α<0.002;

[0017] FIG. 3 is a graph showing the plot of L_σ(x₁,α) for 0<α<0.002;

[0018] FIGS. 4 to 42 are plots of optimal α for dimensions 0-38, respectively, of leaf 513;

[0019] FIG. 43 is a graph showing a comparison of α update versus nonupdate for α=1;

[0020] FIG. 44 is graph showing a comparison of α update versus nonupdate for α=2;

[0021] FIG. 45 is a graph of log-likelihood gains of α update formula as a function of iteration;

[0022] FIG. 46 is block diagram of a computer system of the type on which the subject invention may be implemented; and

[0023] FIGS. 47A and 47B, taken together, are a flow diagram illustrating the logic of the computer implemented process according to the invention.

[0024] For the specific problem of automatic machine recognition of speech as considered by S. Basu and C. Micchelli, supra, preliminary attempts were made to find the “best” value of α by running the speech recognizer for several different values of α and then choosing the value of α corresponding to the lowest word recognition error. The choice of recognition accuracy as the optimality criterion for α is dictated by the practical need for best recognizer. While other criteria for recognition accuracy (e.g., error rates based on phonemes, syllables or at a more detailed level of grannularity, the classification accuracy of feature vectors) can be thought of, we are primarily concerned with automated methods of finding an optimal value of the parameter α.

[0025] In the general context of multivariate mixture densities, this is a difficult task; therefore, the description which follows is restricted to the case of modeling univariate data with one mixture component of this form. In this setting, we present an assortment of iterative formulas for estimating the parameter α from the data, and study their numerical performance. For mixtures of multivariate α-densities, we resort to a technique arising in the context of EM estimation methodology. Recall from S. Basu and C. A. Micchelli, supra, and Christopher M. Bishop, supra, that the increase in the log-likelihood is bounded from below (via Jensen's inequality) by a term (often referred to as the Q-function in the EM literature) that can be conveniently maximized via optimization techniques. In the case of the present invention, we carry out the maximization with respect to the extra parameter α. In this regard, we take advantage of the special dependency of Q on α. Moreover, we extend this method to include modeling data by using mixture components each having different values of α.

[0026] Our experimental results in the context for modeling acoustic feature vectors for the recognition of speech indicate that optimal values of α are smaller than one. In this way, we further substantiate the non-Gaussian form of the mixture components for modeling of speech feature vectors in acoustic space.

[0027] 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).$

[0028] Assume that the data comes from a parametric family p(·|λ), λ∈Ω⊂R^qof probability densities. Maximum likelihood estimation (MLE) method demands that the log-likelihood $\begin{array}{cc}{L}_N\left(\lambda \right)=\sum _{i=1}^N\log p\left(x^i❘\lambda \right)& \left(6\right)\end{array}$

[0029] as a function of λ is made as large as is practical. This specifies a data dependent choice of λ and hence a density function to model the data. When N is large and the data is selected as random samples from the density p(·|θ), the quantity 1/NL_Napproaches

[0030] E_θ(log p(·|λ).  (7)

[0031] By Jensen's inequality $\begin{array}{c}{E}_{\theta }\left\{\log p\left(·❘\lambda \right)\right\}-E_{\theta }\left\{\log p\left(·❘\theta \right)\right\}=E_{\theta }\left\{\log \frac{p\left(·❘\lambda \right)}{p\left(·❘\theta \right)}\right\}\\ \ge \log {\int }_{R^d}\frac{p\left(x❘\lambda \right)}{p\left(x❘\theta \right)}p\left(x❘\theta \right)dx\\ =0\end{array}$

[0032] Therefore, the global maximum of (7) occurs for λ=θ. We add to this fact by observing if equality occurs above, i.e.,

[0033] E_θ{log p(·|λ)}=E_θ{log p(·|θ)},  (8)

[0034] then p(x|λ)=p(x|θ), a.e., x∈R^d(provided every member of the family vanishes on sets of measure zero). Thus, whenever any probability density in our family determines its parameter uniquely, we conclude that λ=θ is the unique global maximum of the function in (7). In particular, the family of elliptical densities of the type (6) have this property. We might add that there are parametric families of importance in applications for which the density does not determine its parameters, for example mixture models of Gaussians.

[0035] In the remainder of this section, we improve upon these facts univariate densities of the form $f\left(t\right)=\frac{\rho }{\sqrt{\sigma }}\exp \left(-h\left(\gamma \frac{{\left(t-\hat{\mu }\right)}^2}{\sigma }\right)\right),t\in R,$

[0036] where $\begin{array}{cc}\gamma =\frac{{\int }_{R_{+}}{t}^{\frac{1}{2}}\exp \left(-h\left(t\right)\right)dt}{{\int }_{R_{+}}{t}^{-\frac{1}{2}}\exp \left(-h\left(t\right)\right)dt}\mathrm{and}& \left(9\right)\\ \rho =\frac{{{\int }_{R_{+}}{t}^{\frac{1}{2}}\exp \left(-h\left(t\right)\right)dt)}^{\frac{1}{2}}}{{{\int }_{R_{+}}{t}^{-\frac{1}{2}}\exp \left(-h\left(t\right)\right)dt)}^{\frac{3}{2}}}.& \left(10\right)\end{array}$

[0037] For this class of densities, our parameter space consists of the mean μ∈R, the variance σ∈R₊ and the function h. We use the notation θ=(μ,σ,h) to denote the parameters that determine this density which we now denote by f(·|θ).

[0038] Lemma 1 Let θ=(μ,σ,h), {circumflex over (θ)}=({circumflex over (μ)},{circumflex over (σ)},ĥ) and suppose h and ĥ are increasing functions of R₊. Then, the function F(μ)=E_{{circumflex over (θ)}}(log f(·|θ)) has a unique local maximum at μ={circumflex over (μ)}.

[0039] Proof: Direct computation confirms the equation. $F\left(\mu \right)=\log \rho -\frac{\log \sigma }{2}-{\int }_{+\infty }^{+\infty }h\left(\frac{\gamma {\left(t-\mu \right)}^2}{\sigma }\right)\frac{\rho }{\sqrt{\hat{\sigma }}}\exp \left(-\hat{h}\left(\frac{\hat{\gamma }{\left(t-\hat{\mu }\right)}^2}{\hat{\sigma }}\right)\right)dt.$

[0040] To prove the claim, we show that F′(μ) is negative for μ>{circumflex over (μ)}. To this end, we observe that $F^{\prime }\left(\mu \right)=\frac{2\gamma \hat{\rho }}{\sigma \sqrt{\hat{\sigma }}}-{\int }_0^{\infty }{\mathrm{th}}^{\prime }\left(\gamma \frac{t^2}{\sigma }\right)\left[g\left({\left(t-\mu +\hat{\mu }\right)}^2\right)-g\left({\left(t-\hat{\mu }+\mu \right)}^2\right)\right]dt,$

[0041] where $g\left(t\right):=\exp \left(-\hat{h}\left(\hat{\gamma }\frac{t}{\hat{\sigma }}\right)\right)t\in R.$

[0042] Let us first point out that F′(μ)<0 for μ>{circumflex over (μ)}. This claim follows by noting that (y−μ+{circumflex over (μ)})²<(y+μ−{circumflex over (μ)})²when μ>{circumflex over (μ)} and y>0 and that g(·) is a decreasing function.

[0043] Lemma 2 Let θ=(μ,σ,h), {circumflex over (θ)}=({circumflex over (μ)},{circumflex over (σ)},ĥ). Suppose that h(x) and xh′(x) are increasing functions of x for x>0. Then the function F(μ,σ)=E_{{circumflex over (θ)}}(log f(·|θ)) has a unique global maximum at μ={circumflex over (μ)}, σ={circumflex over (σ)}.

[0044] Proof: In Lemma 1, we proved that for any σ the function F(·,σ) is maximized for μ={circumflex over (ν)}. Thus, it suffices to consider the function G:=F ({circumflex over (μ)},·). The function has the form $G\left(\sigma \right)=\log \rho -\frac{\log \sigma }{2}-2\rho {\int }_0^{\infty }h\left(\gamma \frac{\hat{\sigma }}{\sigma }{y}^2\right)\exp \left(h\left(\gamma y^2\right)\right)dy,\sigma \in R_{+}$

[0045] Alternatively, we have the equation $G\left(\sigma \right)=\log \rho +\log \phi -\frac{1}{2}\log \hat{\sigma }-2\rho {\int }_0^{\infty }\tilde{h}\left(\phi y\right)\exp \left(-\tilde{h}\left(y\right)\right)dy.$

[0046] where {tilde over (h)}(y):=h(γy²), y∈R₊ and $\phi =\sqrt{\frac{\hat{\sigma }}{\sigma }}.$

[0047] Therefore, we obtain $\begin{array}{cc}{G}^l\left(\sigma \right)=-\frac{x^3}{2\hat{\sigma }}\left\{\frac{1}{x}-2\rho {\int }_0^{\infty }y{\tilde{h}}^l\left(\mathrm{xy}\right)\exp \left(-\tilde{h}\left(y\right)\right)dy\right\}.& \left(11\right)\end{array}$

[0048] Moreover, using integration by parts, we get $\begin{array}{cc}\begin{array}{c}2\rho {\int }_0^{\infty }y{\tilde{h}}^l\left(y\right)\exp \left(-\tilde{h}\left(y\right)\right)dy=2\rho {\int }_0^{\infty }y\frac{d}{dy}\left(\exp \left(-\tilde{h}\left(y\right)\right)\right)dy\\ =2\rho {\int }_0^{\infty }\exp \left(-\tilde{h}\left(y\right)\right)dy\\ =\rho {\int }_{-\infty }^{\infty }\exp \left(-\tilde{h}\left(y\right)\right)dy\\ =1\end{array}& \left(12\right)\end{array}$

[0049] Combining these two equations, we have that $G^l\left(\sigma \right)=-\frac{\rho x^2}{\hat{\sigma }}{\int }_0^{\infty }\exp \left(-\tilde{h}\left(y\right)\right)\left\{y{\tilde{h}}^l\left(y\right)-xy{\tilde{h}}^l\left(xy\right)\right\}.$

[0050] By hypothesis, we conclude that (y{tilde over (h)}^l(y)−xy{tilde over (h)}^l(xy)) is negative if x>1 and, positive if x<1. This implies that G′(σ) is negative if σ>{circumflex over (σ)} and G′(σ) is positive if σ<{circumflex over (σ)} thereby proving the claim.

[0051] The next result is restricted to the family of functions

[0052] h(t)=t^α/2, α>0t∈R₊

[0053] In this case, our parameter vector has the form θ=(μ,σ,α). To state the next result, we observer that the constants in equations (9) and (10) are given by $\begin{array}{cc}\gamma \left(\alpha \right)=\frac{\Gamma \left(3/\alpha \right)}{\Gamma \left(1/\alpha \right)}\mathrm{and}\rho \left(\alpha \right)=\frac{{\mathrm{\alpha \Gamma }}^{1/2}\left(3/\alpha \right)}{2{\Gamma }^{3/2}\left(1/\alpha \right)}.& \left(13\right)\end{array}$

[0054] Lemma 3 Suppose that θ=(μ,σ,α), {circumflex over (θ)}=({circumflex over (μ)},{circumflex over (σ)},{circumflex over (α)}). Let H(μ,σ)=E_{{circumflex over (θ)}}(log f(·|θ)). Then H has a unique global maximum at μ={circumflex over (μ)}, σ={circumflex over (σ)} where $\beta \left(\alpha ,\hat{\alpha }\right)={\left\{\frac{\mathrm{\alpha \Gamma }\left(\frac{\alpha +1}{\hat{\alpha }}\right)}{\Gamma \left(\frac{1}{\hat{\alpha }}\right)}\right\}}^{\frac{2}{\alpha }}\frac{\Gamma \left(\frac{3}{\alpha }\right)\Gamma \left(\frac{1}{\hat{\alpha }}\right)}{\Gamma \left(\frac{3}{\hat{\alpha }}\right)\Gamma \left(\frac{1}{\alpha }\right)}$

[0055] Proof: It follows from Lemma 1 that for any σ>0, the function H(·,σ) is maximized for μ={circumflex over (μ)}. Furthermore, by following the proof of Lemma 2, we get $H\left(\hat{\mu },\sigma \right)=\log \rho \left(\alpha \right)-\frac{1}{2}\log \sigma \left(\alpha \right)-\frac{2\rho \left(\hat{\alpha }\right)}{{\gamma \left(\hat{\alpha }\right)}^{1/2}}{\int }_0^{\infty }{\left(\frac{\gamma \hat{\alpha }}{\gamma \left(\hat{\alpha }\right)\sigma }{t}^2\right)}^{\alpha /2}\exp \left(-t^{\hat{\alpha }}\right)dt.$

[0056] Evaluating the above integral, we obtain the equation $H\left(\hat{\mu },\sigma \right)=\log \rho -\frac{1}{2}\log \sigma -\frac{2\rho \hat{\sigma }}{{\gamma \left(\hat{\alpha }\right)}^{1/2}}{\left(\frac{\gamma \hat{\alpha }}{\gamma \left(\hat{\alpha }\right)\sigma }\right)}^{\alpha /2}\frac{1}{\hat{\alpha }}\Gamma \left(\frac{\alpha +1}{\hat{\alpha }}\right).$

[0057] We proceed by differentiating the function H with respect to σ. $\frac{dH\left(\hat{\mu },\sigma \right)}{d\sigma }=\frac{1}{\sigma }\left\{-\frac{1}{2}+\frac{\alpha }{\hat{\alpha }}\frac{\rho \left(\hat{\alpha }\right)}{{\gamma \left(\hat{\alpha }\right)}^{1/2}}{\left(\frac{\gamma \hat{\sigma }}{\gamma \left(\hat{\alpha }\right)\sigma }\right)}^{\alpha /2}\Gamma \left(\frac{\alpha +1}{\hat{\alpha }}\right)\right\}$

[0058] It is now easily verified that $\frac{dH\left(\hat{\mu },\sigma \right)}{d\sigma }=\{\begin{array}{cc}<0,& \mathrm{if}\sigma >\beta \left(\alpha ,\hat{\alpha }\right)\hat{\sigma }\\ >0,& \mathrm{if}\sigma <\beta \left(\alpha ,\hat{\alpha }\right)\hat{\sigma }\end{array}$

[0059] Thus, σ=ψ(α){circumflex over (σ)} is the global maximum of the function H({circumflex over (μ)},·).

[0060] Remark: Note that we have β(α,{circumflex over (α)}))=1.

[0061] To prove the main result of the section, we need a fact about the gamma function. As we will need to differentiate log Γ(x), we introduce the notation $\Psi \left(x\right)=\frac{d}{dx}\left(\log \Gamma \left(x\right)\right)=\frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}$

[0062] The function Ψ(x) is known as the digamma function and Ψ′(x) as the trigamma function. From M. Abramowitz and I. Stegan, Handbook of Mathematical Statistics, Dover Publications, New York, Ninth Dover printing (1972), the trigamma function has the representation $\begin{array}{cc}{\Psi }^{\prime }\left(x\right)=\sum _{n=0}^{\infty }\frac{1}{{\left(x+n\right)}^2}\mathrm{for}x>0& \left(14\right)\end{array}$

[0063] and satisfies Ψ′(x)>0 for x>0 and $\Psi \left(x+1\right)=\frac{1}{x^2}+\Psi \left(x\right).$

[0064] We shall use this in proving the following theorem.

[0065] Theorem 1 Suppose that θ=(μ,σ,α) and {circumflex over (θ)}=({circumflex over (μ)},{circumflex over (σ)},{circumflex over (α)}). Let H(μ,σ,α)=E_{{circumflex over (θ)}}{f(·|θ)}. Then, H has a unique local maximum for μ={circumflex over (μ)}, σ={circumflex over (σ)} and α={circumflex over (α)}.

[0066] Proof: From Lemma 1 and Lemma 3, it suffices to consider μ={circumflex over (μ)} and α=β(α,{circumflex over (α)}){circumflex over (σ)}.

[0067] We shall write H(α) for the expression H({circumflex over (μ)},β(α,{circumflex over (α)}){circumflex over (σ)},α) and get $H\left(\alpha \right)=\log \rho -\frac{1}{2}\log \beta \left(\alpha ,\hat{\alpha }\right)\hat{\sigma }-\frac{2\hat{\rho }}{\hat{\alpha }\sqrt{\hat{\gamma }}}{\left(\frac{\gamma \hat{\sigma }}{\beta \left(\alpha ,\hat{\alpha }\right)\hat{\gamma }\hat{\sigma }}\right)}^{\alpha /2}\Gamma \left(\frac{\alpha +1}{\hat{\alpha }}\right).$

[0068] We define t=1/α, {circumflex over (t)}=1/{circumflex over (α)}, {overscore (H)}(t)=H(α) and observe that

[0069] {overscore (H)}(t)=−log Γ(1+t)−t+t+log Γ({circumflex over (t)}+{circumflex over (t)}/t)+tlog Γ({circumflex over (t)})+f½ log({circumflex over (γ)}/{circumflex over (σ)})−log 2.

[0070] Our goal is to show that for every {circumflex over (t)}>0 the function {overscore (H)} has a unique local maximum for t>0 at t={circumflex over (t)}. Before doing so, let us observe that $\begin{array}{cc}\underset{t\to \infty }{\lim }\stackrel{\_}{H}\left(t\right)=-\infty & \left(15\right)\end{array}$

[0071] for the proof of this claim, we invoke Stirling's formula (see Abramowitz and Stegan, supra)

[0072] log Γ(t+1)=(t+½)log t−(t+1)+1/2 log 2π+O(1/t).

[0073] Substituting this relation into the definition of {overscore (H)} we get $\begin{array}{c}\stackrel{\_}{H}\left(t\right)=-\left(\left(t+\frac{1}{2}\right)\log t-\left(t+1\right)+\frac{1}{2}\log 2\pi \right)-t+t\log t-\\ \left(t\log \Gamma \left(\hat{t}\right)+t\frac{\hat{t}}{t}\Psi \left(\hat{t}\right)\right)+t\log \Gamma \left(\hat{t}\right)+\frac{1}{2}\log \left(\frac{\hat{\gamma }}{\hat{\sigma }}\right)-\log 2+O\left(\frac{1}{t}\right)\\ =-\frac{1}{2}\log t+1-\frac{1}{2}\log 2\pi -\hat{t}\Psi \left(\hat{t}\right)+\frac{1}{2}\log \left(\frac{\hat{\gamma }}{\hat{\sigma }}\right)-\log 2+O\left(\frac{1}{t}\right).\end{array}$

[0074] from which (15) follows. Returning to our goal to show that {overscore (H)} has a unique local maximum at {overscore (t)}, we compute the derivative ${\stackrel{\_}{H}}^l\left(t\right)=-\Psi \left(t+1\right)+\frac{\hat{t}}{t}\Psi \left(\hat{t}+\frac{\hat{t}}{t}\right)+\log \left(\frac{t}{\hat{t}}\frac{\Gamma \left(\hat{t}+1\right)}{\Gamma \left(\hat{t}+\frac{\hat{t}}{t}\right)}\right)$

[0075] from which we can see that {overscore (H)}^l({circumflex over (t)})=0. The second derivative of {overscore (H)} is given by $\begin{array}{cc}{\stackrel{\_}{H}}^{\mathrm{ll}}\left(t\right)=\frac{1}{t}-{\Psi }^l\left(t+1\right)-\frac{{\hat{t}}^2}{t^3}{\Psi }^l\left(\hat{t}+\frac{\hat{t}}{t}\right).& \left(16\right)\end{array}$

[0076] We shall now demonstrate that the second drivative has a unique simple zero for t>0. The first step in this argument is to show that {overscore (H)}^llis strictly negative for t∈[0,2{circumflex over (t)}]. To this end, we require the following inequality $\begin{array}{cc}{\Psi }^l\left(t\right)>\frac{1}{t}+\frac{1}{2t^2},t>0.& \left(17\right)\end{array}$

[0077] For the proof, we use the convexity of the function t→t⁺²and the trapezoidal rule, which gives the estimate $\begin{array}{cc}\begin{array}{c}\sum _{n=0}^{\infty }\frac{1}{2}\left\{\frac{1}{{\left(t+n\right)}^2}+\frac{1}{{\left(t+n+1\right)}^2}\right\}>\sum _{n=0}^{\infty }{\int }_{t+n}^{t+n+1}\frac{1}{x^2}dx\\ ={\int }_t^{\infty }\frac{1}{x^2}dx\\ =\frac{1}{t}.\end{array}& \left(18\right)\end{array}$

[0078] Substituting equation (14) into the above inequality yields the desired conclusion. Next, we use (17) and the expression (16) for {overscore (H)}¹¹to obtain the inequality ${\stackrel{\_}{H}}^{\mathrm{ll}}\left(t\right)\le \frac{t-2\hat{t}}{2t^2\left(t+1\right)},t>0.$

[0079] This inequality is insufficient to prove the result. Further estimates for HΔ are required.

[0080] The preceding lemma indicates the possibility that μ={circumflex over (μ)}, σ={circumflex over (σ)}. and α={circumflex over (α)} may indeed be the global maximum for H(μ,σ,α). We conjecture that this is the case. However, not knowing the true value of {circumflex over (μ)}, {circumflex over (σ)} and {circumflex over (α)} we can, of course, not compute H(μ,σ,α) and are left with maximizing $L\left(\mu ,\sigma ,\alpha \right)=\prod _{i=1}^{N}f\left(x^i;\mu ,\sigma ,\alpha \right).$

[0081] This is equivalent to maximizing $\frac{1}{N}\log L=\log \rho -\frac{1}{2}\log \sigma -{\left(\frac{\gamma }{\sigma }\right)}^{\alpha /2}\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }.$

[0082] One would expect that 1/N log L(μ,σ,α)≈H(μ,σ,α) and, therefore, that the maximum would occur at μ={circumflex over (μ)}, σ={circumflex over (σ)} and α={circumflex over (α)} for large values of N. However, as we show later, max_μ,σ1/N log L(μ,σ,α) goes to infinity as α→0. Recall that this was not so for H(μ,σ,α). One must therefore take care that the value of α does not become too small when seeking a local maximum of $\frac{1}{N}\log L.$

[0083] Towards showing this behavior of $\frac{1}{N}\log L,$

[0084] we find a lower bound for ${\max }_{\mu ,\sigma }\frac{1}{N}\log L\left(\mu ,\sigma ,\alpha \right).$

[0085] Lemma 4 Let $L_{\sigma }\left(\mu ,\alpha \right)={\max }_{\sigma \ge 0}\frac{1}{N}\log L\left(\mu ,\sigma ,\alpha \right).\mathrm{If}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }>0$

[0086] then $\begin{array}{cc}{L}_{\sigma }\left(\mu ,\alpha \right)=\log \left(\frac{\alpha }{2}\right)-\frac{1}{\alpha }\log \alpha -\log \Gamma \left(\frac{1}{\alpha }\right)-\frac{1}{\alpha }-\frac{1}{\alpha }\log \left(\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)& \left(19\right)\end{array}$

[0087] Proof: Evaluating log L yields the expression $\frac{1}{N}\log L\left(\mu ,\sigma ,\alpha \right)=\log \rho -\frac{1}{2}\log \sigma -{\left(\frac{\gamma }{\alpha }\right)}^{\frac{\alpha }{2}}\frac{1}{2}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }.$

[0088] Clearly, $\frac{1}{N}\log L\to -\infty$

[0089] as σ→0 and σ→∞. Since log L is continuously differentiable with respect to σ for 0<σ<∞, it follows that ${\frac{\partial }{\partial \sigma }\left[\left(\frac{1}{N}\log L\right)\right]}_{\sigma ={\sigma }_0}=0$

[0090] at the maximizing value for σ, namely σ₀. Differentiating, we get: $\frac{\partial }{\partial \sigma }\left(\frac{1}{N}\log L\right)=-\frac{1}{2\sigma }+\frac{\alpha }{2\sigma }{\left(\frac{\gamma }{\sigma }\right)}^{\frac{\alpha }{2}}\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }$

[0091] Solving for $\frac{\partial }{\partial \sigma }\left(\frac{1}{N}\log L\right)=0$

[0092] we get the equations $\begin{array}{cc}{\left(\frac{\gamma }{{\sigma }_0}\right)}^{\frac{\alpha }{2}}\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }=\frac{1}{\alpha }\mathrm{and}& \left(20\right)\\ {\sigma }_0={\gamma \left(\frac{\alpha }{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)}^{\frac{2}{\alpha }}& \left(21\right)\end{array}$

[0093] We now have $L_{\sigma }\left(\mu ,\alpha \right)=\log \rho -\frac{1}{2}\log {\sigma }_0-{\left(\frac{\gamma }{{\sigma }_0}\right)}^{\frac{\alpha }{2}}\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }$

[0094] Substituting in (20) we get that $L_{\sigma }\left(\mu ,\alpha \right)=\log \rho -\frac{1}{2}\log {\sigma }_0-\frac{1}{\alpha }$

[0095] and (21) gives $\begin{array}{cc}{L}_{\sigma }\left(\mu ,\alpha \right)=\log \left(\frac{\rho }{\sqrt{\gamma }}\right)-\frac{1}{\alpha }\log \left(\frac{\alpha }{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)& \left(22\right)\end{array}$

[0096] Equation (19) is then arrived at by using (13) and (22).

[0097] After maximizing over σ, one needs to maximize over μ. The problem with maximizing over μ, however, is that for α<1 the derivative with respect to μ of 1/N log L(μ,σ₀,α) develops singularities at the data points and the maximum may occur at any of the data points as well as at value of μ satisfying $\frac{d}{d\mu }\left\{\log L\left(\mu ,{\sigma }_0,\alpha \right)\right\}=0.$

[0098] We will show that if μ=x^jfor some j∈{1,2, . . . ,N} then log L(x^j,σ₀,α)→∞ as α→0. To achieve this, we must understand what happens to $\frac{1}{\alpha }\log \left(\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)\mathrm{as}\alpha \to 0.$

[0099] Now the arithmetic geometric inequality becomes and equality as α→0, i.e., $\underset{\alpha \to 0}{\lim }{\left(\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)}^{\frac{1}{\alpha }}=\prod _{i=1}^N{\left|x^i-\mu \right|}^{\frac{1}{N}}$

[0100] if |xⁱ−μ|>0 for all i=1,2, . . . ,N (see page 15 of G. Hardy, J. E. Littlewood and G. Polya, Inequalities, Cambridge Mathematical Library, 1991).

[0101] For small values of α, the approximation $\frac{1}{\alpha }\log \left(\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)=\frac{1}{N}\sum _{i=1}^N\log \left|x^i-\mu \right|+O\left(\alpha \right)$

[0102] holds if |xⁱ−μ|>0 for i=1,2, . . . ,N, considering the more general case, where μ=x^jfor some i∈{1,2, . . . , N}, we introduce the set S={i: x_i≠μ_i, i∈{1,2, . . . ,N}}. We then have $\begin{array}{c}\frac{1}{\alpha }\log \left(\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }\right)=\frac{1}{\alpha }\log \left(\frac{1}{N}\sum _{i\in S}{\left|x^i-\mu \right|}^{\alpha }\right)\\ =\frac{1}{\alpha }\log \left(\frac{\left|S\right|}{N}\right)+\left(\frac{1}{\left|S\right|}\sum _{i\in S}\log {\left|x^i-\mu \right|}^{\alpha }\right)+O\left(\alpha \right)\end{array}$

[0103] Using this in addition to Lemma 4 we arrive at Lemma 5.

[0104] Lemma 5 For small values of α we have $L_{\sigma }\left(\mu ,\alpha \right)\ge \frac{1}{2}\log \alpha -\frac{1}{\alpha }\log \left(\frac{\left|S\right|}{N}\right)+\frac{1}{\left|\alpha \right|}\sum _{i\in S}\log \left|x^i-\mu \right|\frac{1}{2}\log \left(\frac{\pi }{2}\right)+O\left(\alpha \right)$

[0105] where S={i: x_i≠μ_i, i∈{1,2, . . . ,N}}.

[0106] Proof: Stirling's formula tells us that $\log \Gamma \left(z\right)=\left(z-\frac{1}{2}\right)\log \left(z\right)-z-\frac{1}{2}\log \left(2\pi \right)+O\left(\frac{1}{z}\right)$

[0107] for large values of z. If α is small, we, therefore, have $\log \Gamma \left(\frac{1}{\alpha }\right)=\left(\frac{1}{\alpha }-\frac{1}{2}\right)\log \left(\frac{1}{\alpha }\right)-\frac{1}{\alpha }-\frac{1}{2}\log \left(2\pi \right)+O\left(\alpha \right)$

[0108] Together with the arithmetic-geometric connection this provides us $\begin{array}{c}{L}_{\sigma }\left(\mu ,\alpha \right)=\log \left(\frac{\alpha }{2}\right)-\frac{1}{\alpha }\log \left(2\pi \right)-\frac{1}{\alpha }-\\ \left\{\left(\frac{1}{\alpha }-\frac{1}{2}\right)\log \alpha -\frac{1}{\alpha }-\frac{1}{2}\log \left(2\pi \right)+O\left(\alpha \right)\right\}-\\ \left\{\frac{1}{\alpha }\log \left(\frac{\left|S\right|}{N}\right)+\frac{1}{\left|S\right|}\sum _{i\in S}^{}\log \left|x^i-\mu \right|+O\left(\alpha \right)\right\}\\ =\frac{1}{2}\log \alpha -\frac{1}{\alpha }\left(\frac{\left|S\right|}{N}\right)+\frac{1}{\left|S\right|}\sum _{i\in S}^{}\log \left|x^i-\mu \right|+\\ \frac{1}{2}\log \left(\frac{\pi }{2}\right)+O\left(\alpha \right).\end{array}$

[0109] Clearly, if μ=xⁱfor some i∈{1,2,3, . . . ,N} we have |S|<N. Thus, fixing μ=x_jwe get $L_{\sigma }\left(\mu ,\alpha \right)\ge \frac{1}{2}\log \alpha -\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right)+\frac{1}{\left|S\right|}\sum _{i\in S}^{}\log \left|x^i-\mu \right|+\frac{1}{2}\log \left(\frac{\pi }{2}\right)+O\left(\alpha \right)$

[0110] But $\frac{1}{\left|S\right|}\sum _{i\in S}\log \left|x^i-\mu \right|$

[0111] is independent of α and log α is dominated by $\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right)$

[0112] and so we see that

[0113] L_σ(μ,α)→∞ for α→∞.

[0114] This was not the case for H(μ,σ,α). On the other hand we see that $-\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right)\text{ starts to dominate}\frac{1}{N}L\left(\mu ,\sigma ,\alpha \right)\text{ when}$ $\frac{1}{2}\log \alpha -\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right)>0$

[0115] We would like a rule of thumb that ensures that $\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right)$

[0116] does not dominate. We pick $\alpha =\frac{1}{N}$

[0117] and evaluate $\frac{1}{2}\log \alpha -\frac{1}{\alpha }\log \left(1-\frac{1}{N}\right){❘}_{\alpha =\frac{1}{N}}=-\frac{1}{2}\log N+1+O\left(\frac{1}{N}\right)$

[0118] But −½ log N+1<0 if N>exp(2)=7.3891, so we can safely use the rule of thumb that α>1/N for our experience as we will mostly that N≧8.

[0119] We conclude this section by noting that the “optimal” choice of α is a local maximum of 1/N log(μ,σ,α) for which α>1/N.

[0120] To verify the theory of preceding sections, we generated synthetic data from a standard Gaussian density, x_i, i=1,2, . . . ,N. Lemma 5 provides an explicit formula for L_σ(μ,α). Plotting this function for μ and σ for our synthetic data, we expect to see a local maximum at μ=0, σ=2 according to Theorem 1. Also, we expect the log-likelihood to go to infinity for small values of α. We make plots of L_σ(μ,α) for values of μ and σ around 0 and 2 respectively, and for small values of α (with α<1/N).

[0121] FIG. 1 and FIG. 2 illustrate these two points. Investigating the value of average log-likelihood in the figures, we find that the maximum likelihood in FIG. 1 is approximately −1.3629 and is attained for μ=−0.06 and α=2.109. In FIG. 2 the likelihood is −3.9758. This is not according to our intention. The reason for this is that FIG. 2 is generated on a rectangular grid, and the data points happen not to fall on the grid points. If we plot L_σ(x_l,α) the situation is different. FIG. 3 shows L_σ(x_l,α) as a function of α. We see that the likelihood is larger than −3.9758 and −1.3629 and indeed seems to grow towards infinity as α goes to zero. This graph confirmed the behavior predicted in Lemma 5. After our theoretical exposition, it seems appropriate to shift gear and consider some real data. We are, in particular, interested in acoustic data of speech. Digitized speech sampled at a rate of 16 KHz is considered. A frame consists of a segment of speech of duration 25 msec. and produces thirty-nine dimensional acoustic cepstral vectors via the following process, which is standard in speech recognition literature. Frames are advanced every 10 msec. to obtain succeeding acoustic vectors.

[0122] First, magnitudes of discrete Fourier transform of samples of speech data in a frame are considered in a logarithmically warped frequency scale. Next, these amplitude values themselves are transformed to a logarithmic scale, and subsequently, a rotation in the form of discrete cosine transform is applied. The first thirteen components of the resulting vector are retained. First and the second order differences of the sequence of vectors so obtained are then appended to the original vector to obtain the thirty-nine dimensional cesptral acoustic vectors.

[0123] As in supervised learning tasks, we assume that these vectors are labeled according to the basic sounds they correspond to. In fact, the set of 46 phonemes are subdivided into a set of one hundred twenty-six different variants each corresponding to a “state” in the hidden Markov model used for recognition proposed. They are further subdivided into more elemental sounds called allophones by using the method of decision trees depending on the context in which they occur. See Frederick Jelenik, “Statistical Methods for Speech Recognition”, MIT Press, 1997, L. R. Bahl, P. V. Desouza, P. S. Gopaikrishnan, M. A. Picheny, Context dependent vector quantization for continuous speech recognition, Proceedings of IEEE Int. Conf. On Acoustics Speech and Signal Processing, pp. 632-635, 1993, and Leo Breiman, “Classification and Regression Trees” Wadsworth International, Belmont, Calif., 1983.

[0124] We chose an allophone at random and made plots of L_σ(μ,α) for each dimension to determine the best value of α. The particular allophone we chose corresponded to the phone AX as in the world “MAXIMUM” pronounced M AB K S A MUM. These plots are shown in FIGS. 4 to 42. Looking at the figures, we observe that L_σ(μ,α) is highly unsymmetric. Moreover, it appears that the first and second difference (dimensions 12 to 25, FIGS. 12 to 29, and dimensions 26 to 38, FIGS. 30 to 42) generally prefer lower values of α than the unappended cepstral vector.

[0125] Table 1 shows the optimal values of α for each dimension corresponding to the phonemes AX, F, N and T. These phonemes were chosen so as to represent four distinct sounds corresponding to vowel, fricative, nasal and stop. The optimal values of α were extracted from two dimensional surface plots of L_σ(μ,α) for the respective allophones and dimensions.

[0126] As seen from the foregoing discussion, the optimal choice for σ is given by (21) for a fixed μ and α. Once this choice of σ is made the average log-likelihood is giving by (19). Therefore, we need only worry about finding μand α. Since we are unable to do this analytically, we aim at finding a numerically convergent iteration scheme that converges to the optimal values. For the value of μ, we use the same approach as in Basu and Micchelli, supra,. Differentiating L_σ(μ,α) with respect to μ we get $\frac{\partial L_{\sigma }\left(\mu ,\alpha \right)}{\partial \mu }=\frac{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha -2}\left(x_i-\mu \right)}{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\alpha }}$

[0127] Equating this to zero provides the stationary equation for μ, which we can rewrite as $\begin{array}{cc}u=\frac{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\left(\alpha -2\right)}{x}_i}{\sum _{i=1}^N{\left|x_i-\mu \right|}^{\left(\alpha -2\right)}}& \left(23\right)\end{array}$

[0128] Considering μ's on the right hand side as old values, we compute updated values of μ from (23). Despite the analysis performed in Basu and Micchelli, supra, the iterative formula for computing μ so obtained, is seen to converge numerically.

[0129] It remains to find the optimal choice of α. To this end, we consider a portfolio of iterative methods, each of which are tested numerically. Clearly, we know a priori that α>0. Any iteration method considered ought to guarantee that this is so regardless of the starting value of α. Differentiating L_σ(μ,α) with respect to α we obtain the stationary equation for α. $\begin{array}{cc}\log \left(\alpha \right)+\psi \left(\frac{1}{\alpha }+1\right)+\log \left(S\left(\alpha \right)\right)-\frac{\alpha S^l\left(\alpha \right)}{S\left(\alpha \right)}=0& \left(24\right)\end{array}$

[0130] where $S\left(\alpha \right)=\frac{1}{N}\sum _{i=1}^N{\left|x^i-\mu \right|}^{\alpha }.$

[0131] Here, we have used the identity $\psi \left(\frac{1}{\alpha }+1\right)=\psi \left(\frac{1}{\alpha }\right)+\alpha .$

[0132] One way of ensuring that α is positive is by isolating log α on one side and exponentiating both sides to get $\begin{array}{cc}\alpha =\frac{\exp \left(\frac{\alpha S^l\left(\alpha \right)}{S\left(\alpha \right)}-\psi \left(\frac{1}{\alpha }+1\right)\right)}{S\left(\alpha \right)}.& \left(25\right)\end{array}$

[0133] Rewriting this gives the alternate iteration $\begin{array}{cc}\alpha ={\alpha }^{2}S\left(\alpha \right)\exp \left(-\frac{\alpha S^l\left(\alpha \right)}{s\alpha }+\psi \left(\frac{1}{\alpha }+1\right)\right).& \left(26\right)\end{array}$

[0134] Using the formula $\psi \left(\frac{1}{\alpha }+1\right)=\psi \left(\frac{1}{\alpha }\right)+\alpha$

[0135] we isolate α as follows $\alpha =\frac{\alpha S^l\left(\alpha \right)}{S\left(\alpha \right)}-\log S\left(\alpha \right)-\log \left(\alpha \right)-\psi \left(\frac{1}{\alpha }\right)$

[0136] The right hand side of the above equation is not necessarily positive. To make it possible we square both sides and divide by α. This gives the iteration formula for α $\begin{array}{cc}\alpha =\frac{1}{\alpha }{\left(\frac{\alpha S^l\left(\alpha \right)}{S\left(\alpha \right)}-\log S\left(\alpha \right)-\log \left(\alpha \right)-\psi \left(\frac{1}{\alpha }\right)\right)}^2.& \left(27\right)\end{array}$

[0137] Modifying this equation, we also get another iteration formula, namely, $\begin{array}{cc}\alpha ={{\alpha }^3\left(\frac{\alpha S^l\left(\alpha \right)}{S\left(\alpha \right)}-\log S\left(\alpha \right)-\log \left(\alpha \right)-\psi \left(\frac{1}{\alpha }\right)\right)}^{-2}.& \left(28\right)\end{array}$

[0138] Finally, we wish to test Newton's method for locating zeros of an equation. However, Newton's method does not guarantee positivity. To circumvent this problem we introduce the temporary variable α=β². We are solving the equation f(α)=0, where $f\left(\alpha \right)=\psi \left(1+\frac{1}{\alpha }\right)+\log \left(\alpha \right)+\log S\left(\alpha \right)-\frac{\alpha S^{\prime }\left(\alpha \right)}{S\left(\alpha \right)}$

[0139] Classical Newton's iteration is $\hat{\alpha }=\alpha -\frac{f\left(\alpha \right)}{f^{\prime }\left(\alpha \right)},$

[0140] where {circumflex over (α)} is the updated estimate of the root. As we introduce α=β²we consider the function g(β)=f(β²). Newton's iteration applied to g(β) gives us $\beta =\beta -\frac{g\left(\beta \right)}{g^{\prime }\left(\beta \right)}.$

[0141] Converting this to a formula involving α we get $\begin{array}{cc}\hat{\alpha }={\left(\sqrt{\alpha }-\frac{f\left(\alpha \right)}{2\sqrt{\alpha }{f}^{\prime }\left(\alpha \right)}\right)}^2& \left(29\right)\end{array}$