MLE–EQUIVARIANCE, DATA TRANSFORMATIONS AND INVARIANT TESTS OF FIT ***Proposed
running head: MLE–equivariance, data transformations, and invariant tests

Muneya Matsui(1)††† Muneya Matsui’s research is partly supported by the JSPS
Grant-in-Aid for Scientific Research C (23K11019)., Simos Meintanis(2),(3)

(1) Department of Business Administration, Nanzan University, Nagoya,
Japan,‡‡‡Email-address<EMAIL_ADDRESS>
(2) Department of Economics, National and Kapodistrian University of Athens,
Athens, Greece§§§On sabbatical leave from the University of Athens, Email-
address<EMAIL_ADDRESS>
(3) Pure and Applied Analytics, North–West University , Potchefstroom, South
Africa

Abstract. We define data transformations that leave certain classes of
distributions invariant, while acting in a specific manner upon the parameters
of the said distributions. It is shown that under such transformations the
maximum likelihood estimators behave in exactly the same way as the parameters
being estimated. As a consequence goodness–of–fit tests based on standardized
data obtained through the inverse of this invariant data–transformation reduce
to the case of testing a standard member of the family with fixed parameter
values. While presenting our results, we also provide a selective review of
the subject of equivariant estimators always in connection to invariant
goodness–of–fit tests. A small Monte Carlo study is presented for the special
case of testing for the Weibull distribution, along with real–data
illustrations.

Keywords. MLE–equivariance; Goodness–of–fit test; Weibull distribution; Test
for exponentiality

##  1 Introduction

Suppose that $\widehat{\theta}_{n}:=\widehat{\theta}_{n}(X_{1},...,X_{n})$ is
the maximum likelihood estimator (MLE) of a parameter (vector) $\theta$ based
on independent copies $(X_{j},\ j=1,...,n)$ of a random variable $X$ with
distribution $F$. Then by the invariance property of the MLE,
$h(\widehat{\theta}_{n})$ is the MLE of the parameter $h(\theta)$ induced by
an one–to–one transformation $h(\cdot)$; see for instance Zacks (1971)[Theorem
5.1.1] or Lin’kov (2005)[Theorem 5.2.1]. (The restriction to one–to–one
transformations is not necessary.) The question is whether
$\vartheta:=h(\theta)$ is the original parameter of another distribution
within the same distributional class and how can we get to this distribution
by a variable transformation $g(X)$ on the random variable $X$, so that the
MLE $\widehat{\vartheta}_{n}$ of $\vartheta$ can be computed, not by employing
the likelihood of the new distribution but by direct reference to the MLE
$h(\widehat{\theta}_{n})$. If this is so we will see that certain
methodological implications emerge that considerably simplify associated
goodness–of–fit (GOF) procedures.

There exist already a couple of well known examples. Let us start with scale
families of distributions, and in this regard assume that $X$ has a
distribution function (DF) $F(x):=F_{c}(x)$ such that $F_{c}(x)=F_{0}(x/c)$,
for some $c>0$, and for some fixed DF $F_{0}$. In this case the
family–preserving variable transformation is given by $Y=g(X)$ with $g(X)=aX$,
and the new parameter corresponding to the transformed random variable $Y$ is
$\sigma=h(c)$ with $h(c)=ac$, for each $a>0$. Then it may be shown that the
MLE satisfies
$\widehat{\sigma}_{n}(Y_{1},...,Y_{n})=a\widehat{c}_{n}(X_{1},...,X_{n})$,
where $Y_{j}=aX_{j},\ (j=1,...,n)$; see e.g. Problem 3.1 (b), p. 211, in
Lehmann and Casella (1998). The other well known example emerges in
location–scale families whereby the DF of $X$ is given by
$F_{\theta}(x)=F_{0}((x-\delta)/c))$, for some $\theta=(\delta,c)$,
($\delta\in\mathbb{R}$ and $c>0$), with family–preserving transformation
$g(X)=aX+b$, ($a>0$ and $b\in\mathbb{R}$). Then the new parameter is
$\vartheta=(\mu,\sigma)$, where $\vartheta=h(\delta,c)$ with
$h(\delta,c)=(a\delta+b,ac)$, and the MLE satisfies
$\widehat{\vartheta}_{n}(Y_{1},...,Y_{n})=(a\widehat{\delta}_{n}+b,a\widehat{c}_{n})$,
where $Y_{j}=aX_{j}+b,\ (j=1,...,n)$. Such estimators that mimic the behavior
of corresponding parameters under specific data–transformations are labeled as
equivariant estimators; see Definition 2.5 in Chapter 3 of Lehmann and Casella
(1998) for a rigorous definition.

Here we go beyond these clearly linear cases and pinpoint reasonably rich
families of distributions within which, specific distributions differ in more
essential ways than mere location and/or scale. Specifically assume that the
distribution of $X>0$, depends on a scale parameter $c>0$ and a shape
parameter $\kappa>0$, and that the corresponding DF is such that

(1.1) $\displaystyle
F_{c,\kappa}(x)=F_{0}\left(\left(\frac{x}{c}\right)^{\kappa}\right),$

for some fixed DF $F_{0}(\cdot)$. Specific families of distributions arise by
different “kernel” DF choices $F_{0}(\cdot)$ and include:

  1. 1.

The Weibull distribution with kernel DF, $F_{0}(x)=1-{\rm{e}}^{-x},\ x>0$.

  2. 2.

The Pareto type I distribution with kernel DF, $F_{0}(x)=1-x^{-1},\ x>1$.

  3. 3.

The Frechet distribution with kernel DF, $F_{0}(x)={\rm{e}}^{-1/x},\ x>0$.

In this paper we identify a root–type transformation as a family–preserving
variable transformation $g(\cdot)$ within families satisfying (1.1), as well
as the corresponding parameter transformation $h(\cdot)$ and show that the MLE
satisfies

(1.2) $\displaystyle\widehat{\vartheta}_{n}=h(\widehat{\theta}_{n}),$

for the specific choice of this transformation. This is done in Section 2.
Then in Section 3 we investigate the impact that (1.2) has on certain
goodness–of–fit tests (GOF) for the families of distributions that satisfy
(1.1). At the same time we also provide a non–technical and selective review
of the related literature, including cases of test invariance for linear
family–preserving variable transformations. The article concludes with a small
simulation study and real–data examples in the case of GOF testing for the
Weibull distribution in Section 4.

##  2 Distribution–preserving data transformations and equivariant estimators

In view of (1.1), notice that

(2.1) $\displaystyle X\sim F_{c,\kappa}\implies aX^{1/b}\sim
F_{ac^{1/b},b\kappa},$

for each $a,b>0$. Hence the one–to–one root–type transformation $X\mapsto
aX^{1/b}$ is a family–preserving variable transformation for each specific
family $F_{0}$ satisfying (1.1), and maps $(c,\kappa)$ to
$(ac^{1/b},b\kappa)$. The following proposition shows that (1.2) holds for the
MLE in such families of distributions, with variable transformation
$g(X)=aX^{1/b}$ and parameter transformation $h(c,\kappa)=(ac^{1/b},b\kappa)$.

###### Proposition 2.1

Let $(X_{j},\ j=1,...,n)$ be independent copies of $X\sim F_{c,\kappa}$, and
assume that the DF of $X$ satisfies (1.1). Assume further that the density
$f_{c,\kappa}$ corresponding to $F_{c,\kappa}$ exists, and also that the MLE
$(\widehat{c}_{n},\widehat{\kappa}_{n})$ of the parameter $(c,\kappa)$ exists.
Then the MLE corresponding to $aX^{1/b}\sim F_{ac^{1/b},b\kappa}$ with the
sample $(aX^{1/b}_{j},\ j=1,...,n)$ is given by
$(a\widehat{c}^{1/b}_{n},b\widehat{\kappa}_{n})$.

Proof. We first calculate from (1.1) the density corresponding to
$F_{c,\kappa}$ as

(2.2) $\displaystyle
f_{c,\kappa}(x)=\frac{{\rm{d}}F_{c,\kappa}(x)}{{\rm{d}}x}=\frac{{\rm{d}}F_{0}((\frac{x}{c})^{\kappa})}{{\rm{d}}x}=\frac{\kappa}{c}\left(\frac{x}{c}\right)^{\kappa-1}f_{0}\left(\left(\frac{x}{c}\right)^{\kappa}\right),$

where $f_{0}(x):=f_{1,1}(x)=\frac{{\rm{d}}F_{0}(x)}{{\rm{d}}x}$.

Then by straightforward calculations, it follows that the likelihood function
corresponding to $X$ with the sample $(X_{j},\ j=1,\ldots,n)$ is given by

(2.3)
$\displaystyle{\cal{L}}(X_{1},\ldots,X_{n};c,\kappa)=\frac{\kappa^{n}}{\prod_{j=1}^{n}X_{j}}\prod_{j=1}^{n}\left(\frac{X_{j}}{c}\right)^{\kappa}f_{0}\left(\left(\frac{X_{j}}{c}\right)^{\kappa}\right).$

Recall now that $Y:=aX^{1/b}\sim F_{ac^{1/b},b\kappa}$, and denote by
$(\widehat{d}_{n},\widehat{\eta}_{n})$ the MLE for $Y$ with the sample
$(Y_{j}=aX^{1/b}_{j},\ j=1,\ldots,n)$ , i.e.

$\max_{ac^{1/b},b\kappa}{\cal{L}}(Y_{1},\ldots,Y_{n};ac^{1/b},b\kappa)={\cal{L}}(Y_{1},\ldots,Y_{n};\widehat{d}_{n},\widehat{\eta}_{n}).$

Then, by the definition of $(\widehat{c}_{n},\widehat{\kappa}_{n})$ and
$(\widehat{d}_{n},\widehat{\eta}_{n})$, we observe that for any $(c,\kappa)$,
the following hold true:

$\displaystyle\prod_{j=1}^{n}\left(\frac{bX^{1-1/b}_{j}}{a}\right)\cdot{\cal{L}}(X_{1},\ldots,X_{n};c,\kappa)$
$\displaystyle=\frac{(b\kappa)^{n}}{{\prod_{j=1}^{n}Y_{j}}}\prod_{j=1}^{n}\left(\frac{Y_{j}}{ac^{1/b}}\right)^{b\kappa}f_{0}\left(\left(\frac{Y_{j}}{ac^{1/b}}\right)^{b\kappa}\right)$
$\displaystyle={\cal{L}}(Y_{1},\ldots,Y_{n};ac^{1/b},b\kappa)$
$\displaystyle\leq{\cal{L}}(Y_{1},\ldots,Y_{n};\widehat{d}_{n},\widehat{\eta}_{n})$
$\displaystyle=\prod_{j=1}^{n}\left(\frac{bX^{1-1/b}_{j}}{a}\right)\frac{(\widehat{\eta}_{n}/b)^{n}}{\prod_{j=1}^{n}X_{j}}\prod_{j=1}^{n}\left(\frac{X_{j}}{(\widehat{d}_{n}/a)^{b}}\right)^{\widehat{\eta}_{n}/b}f_{0}\left(\left(\frac{X_{j}}{(\widehat{d}_{n}/a)^{b}}\right)^{\widehat{\eta}_{n}/b}\right)$
$\displaystyle=\prod_{j=1}^{n}\left(\frac{bX^{1-1/b}_{j}}{a}\right)\cdot{\cal{L}}(X_{1},\ldots,X_{n};(\widehat{d}_{n}/a)^{b},\widehat{\eta}_{n}/b)$
$\displaystyle\leq\prod_{j=1}^{n}\left(\frac{bX^{1-1/b}_{j}}{a}\right)\cdot{\cal{L}}(X_{1},\ldots,X_{n};\widehat{c}_{n},\widehat{\kappa}_{n}).$

Now putting $(c,\kappa)=(\widehat{c}_{n},\widehat{\kappa}_{n})$ in the first
line of the argument above we obtain

${\cal{L}}(X_{1},\ldots,X_{n};\widehat{c}_{n},\widehat{\kappa}_{n})={\cal{L}}(X_{1},\ldots,X_{n};(\widehat{d}_{n}/a)^{b},\widehat{\eta}_{n}/b).$

Hence,
$(\widehat{d}_{n},\widehat{\eta}_{n})=(a\widehat{c}_{n}^{1/b},b\widehat{\kappa}_{n})$
should be concluded.

It should be pointed out that there exist other data–transformations that lead
to parameter–free tests, such as the log–linear transformation $\kappa(\log
X-\log c)$ that turns a Weibull variate to a variate following a standard
extreme–value distribution. The difference with the root–type transformation
suggested herein is that this transformation is not distribution–specific, but
instead it applies to rich families of distributions. In this connection, more
general distributions with extra parameters may be included in our framework
of families of distributions. For instance the Burr type XII distribution
belongs to the class of distributions defined by (1.1), with DF
$F_{0}(x)=1-(1+x)^{-\xi},\ \xi>0$, and thus also satisfies Prop. 2.1.
Specifically if $X$ follows a Burr type XII distribution with parameters
$(c,\kappa,\xi)$ then $aX^{1/b}$ follows a Burr type XII distribution with
parameters $(ac^{1/b},b\kappa,\xi)$. Two other well known classes are the
exponentiated Weibull distribution with $F_{0}(x)=(1-{\rm{e}}^{-x})^{\xi},\
\xi>0$, and the generalized gamma distribution with density given by (2.2)
where $f_{0}(x)=(\Gamma(\xi))^{-1}x^{\xi-1}{\rm{e}}^{-x},\ \xi>0$.

To the best of our knowledge the most general class of distributions that
satisfy (1.1) and Prop. 2.1 is the “interpolating family” of distributions on
$(0,\infty)$ recently introduced by Sinner et al. (2023). As it will be seen
below Prop. 2.1 allow us to carry out a much simpler GOF test that refers to a
subfamily of the family under test whereby two (out of two, three or even
four) parameters have been removed.

##  3 Parameter–free test procedures

It should be noted that the notion of equivariant estimators, and invariant
tests (to be discussed below) is a recurring theme in Statistics, with Chapter
3 (resp. Chapter 6) of the classical treatment of Lehmann and Casella (1998)
(resp. Lehmann and Romano (2005)) devoted to such estimators (resp. tests).
Nevertheless the approach in Lehmann and Casella (1998) is mostly based on
estimation optimality, which is not always relevant in the context of GOF
testing. By way of example, smooth tests of fit originally introduced by
Neyman (1937), and more recently studied by Ledwina (1994) and in the
monograph by Rayner et al. (2009), are intrinsically related to moment
estimators which are more often than not less efficient than other estimators.
In fact Klar (2000) points out that the method of moments is the only
meaningful estimation method in the context of smooth tests of fit. Another
point in case is made by Drost et al. (1990) who argue that rather than
estimation optimality, robustness or more precisely luck of it, is important
in the context of GOF testing. Moreover, the test optimality approach often
adopted by Lehmann and Romano (2005)[Chapter 6] (see also Vexler and Hutson
(2023)) is not even feasible in any reasonably wide context of testing, such
as GOF testing with unspecified parameters. This is noted in Lehmann and
Romano (2005)[§14.1], and is formally stated and shown by Janssen (2000) and
Escanciano (2009). (It should be pointed out however that within the narrow
context of testing a distribution against a specific alternative, likelihood
ratio tests applied on maximal invariants lead to optimal invariant tests).
Therefore, it appears that the methodological implications of equivariant
estimators on GOF testing with estimated parameters have not been put forward
beyond the simple linear transformation case of scale or location–scale
families, and even in those cases they have not been sufficiently emphasized.

On the basis of the preceding discussion we motivate our current
parameter–free procedures by starting again with simple scale families of
distributions with DF $F(x)=F_{0}(x/c)$. As already mentioned in the
Introduction, the MLE of the scale parameter $c$ satisfies
$\widehat{c}_{n}(aX_{1},...,aX_{n})=a\widehat{c}_{n}(X_{1},...,X_{n})$. As a
result, any GOF test for such families of distributions, say
$T_{n}(X_{1},...,X_{n})$, that depends on $X_{j}$ only through
$\widehat{Y}_{j}=X_{j}/\widehat{c}_{n}$, $(j=1,...,n)$, is scale invariant,
i.e. it satisfies

(3.1) $\displaystyle T_{n}(aX_{1},...,aX_{n})=T_{n}(X_{1},...,X_{n}),$

for each $a>0$, and consequently the test $T_{n}$ may be applied by assuming
without loss of generality that we are testing for $F=F_{0}$ with $c=1$. (If
$X>0$, this is a special case of Prop. 2.1 for $\kappa=1$, and eqn. (3.1)
follows from eqn. (3.6) below for $b=1$). The by far most popular such testing
problem is that of testing for exponentiality, and the reader is referred to
the review articles of Henze and Meintanis (2005) and Allison et al. (2017),
for scale invariant tests for exponentiality.

One level up are location–scale families whereby the corresponding
location–scale invariant test results by considering the MLE–standardized
observations $\widehat{Y}_{j}=(X_{j}-\widehat{\delta}_{n})/\widehat{c}_{n},\
(j=1,...,n)$, and analogously satisfies

(3.2) $\displaystyle T_{n}(aX_{1}+b,...,aX_{n}+b)=T_{n}(X_{1},...,X_{n}),$

for each $a>0$ and $b\in\mathbb{R}$, implying that $T_{n}$ can be performed by
setting $(\delta,c)=(0,1)$. The monograph of Thode (2002) focuses on the most
popular testing problem here, i.e. that of testing for normality, but without
particular reference to test invariance. (Nevertheless most normality tests
are routinely applied on $(X_{j}-\widehat{\delta}_{n})/\widehat{c}_{n},\
(j=1,...,n)$, with $\widehat{\delta}_{n}$ (resp. $\widehat{c}_{n}$) being the
sample mean (resp. sample standard deviation), which automatically implies the
location–scale invariance stated in (3.2)). On the other hand, Epps (2005)
considers GOF tests for general location–scale families with implicit
reference to invariance. It should be mentioned that the methodology in Epps
(2005) is confined to tests utilizing the empirical characteristic function as
their main tool, which might seem as a somewhat less well known approach, but
with minor modifications the location–scale invariance argued in Epps (2005)
applies more generally to any given GOF test. Moreover, the paper itself as
well as earlier (see e.g., Epps and Pulley (1983), Epps (1993)) and subsequent
works (see e.g., Hall et al. (2013)), including some of the papers to be
referenced herein (see e.g., Meintanis and Swanepoel (2007), Meintanis et al.
(2015)) make a good case about using the empirical characteristic function for
GOF testing, instead of more standard tools such as the empirical DF.

We will briefly digress from univariate distributions, to discuss the very
important case of multivariate data. In this connection we note that
location–scale equivariance has been extended to vectorial observations as
“affine–equivariance”. The most relevant context for affine equivariant
estimators and affine invariant tests is that of (multivariate) elliptical
distributions and an excellent discussion of such estimators and tests may be
found in Hallin and Jurečková (2012). In fact Hallin and Jurečková (2012)
argue that affine–invariant tests, i.e. tests that for arbitrary dimension
$p\geq 1$ satisfy, $T_{n}(AX_{1}+b,...,AX_{n}+b)=T_{n}(X_{1},...,X_{n})$, for
any non–singular $p\times p$ matrix $A$ and any vector $b\in\mathbb{R}^{p}$,
should be based on the Mahalanobis distances between the observed
$p$–dimensional vectors $(X_{j},\ j=1,...,n)$, a point also made by Henze
(2002) in the context of testing for multivariate normality, and by Meintanis
et al. (2015) for the more general elliptically symmetric stable–Paretian
distribution.

Let us return now to our main problem. In this regard, Prop. 2.1 entails that
the MLE estimator $(\widehat{c}_{n},\widehat{\kappa}_{n})$ of $(c,\kappa)$
mimics the equivariance properties of the respective parameters shown in
(2.1), i.e.

(3.3)
$\displaystyle\widehat{c}_{n}(aX^{1/b}_{1},...,aX^{1/b}_{n})=a\>\widehat{c}^{1/b}_{n}(X_{1},...,X_{n}),$

and

(3.4)
$\displaystyle\widehat{\kappa}_{n}(aX^{1/b}_{1},...,aX^{1/b}_{n})=b\>\widehat{\kappa}_{n}(X_{1},...,X_{n}),$

for each $a,b>0$. As a result, any GOF procedure that depends on $(X_{j},\
j=1,...,n)$ only via

(3.5)
$\displaystyle\widehat{Y}_{j}=\left(\frac{X_{j}}{\widehat{c}_{n}}\right)^{\widehat{\kappa}_{n}},\
j=1,...,n,$

satisfies

(3.6) $\displaystyle
T_{n}(aX^{1/b}_{1},...,aX^{1/b}_{n})=T_{n}(X_{1},...,X_{n}),$

for each $a,b>0$, and consequently and without loss of generality, we may
perform the test by assuming $c=\kappa=1$. Clearly, in view of the
distributional invariance of $aX^{1/b}$ figuring in (2.1), test invariance as
illustrated by (3.6) feels like a desirable, even natural, property within the
families of distributions satisfying (1.1). On the practical level eqn. (3.6)
implies that a potentially much simpler test may be invoked for testing
families satisfying (1.1), such as in the case of the Weibull distribution
where any test for exponentiality applied on $(\widehat{Y}_{j},\ j=1,...,n)$
can be used.

###### Remark 3.1

Notice that the data–transformation figuring in eqn. (3.5) is the inverse of
the distribution–preserving root–type transformation $g(X)$ shown in (2.1),
and thus not–surprisingly, its application has a stabilizing effect on the
estimators. Specifically, by replacing in (3.3)–(3.4), $(a,b)$ by
$(\widehat{c}_{n}^{-\widehat{\kappa}_{n}},\widehat{\kappa}^{-1}_{n})$ we see
easily that

$\widehat{c}_{n}(\widehat{Y}_{1},...,\widehat{Y}_{n})=\widehat{c}_{n}\left(\left(\frac{X_{1}}{\widehat{c}_{n}}\right)^{\widehat{\kappa}_{n}},...,\left(\frac{X_{n}}{\widehat{c}_{n}}\right)^{\widehat{\kappa}_{n}}\right)=1,$

and

$\widehat{\kappa}_{n}(\widehat{Y}_{1},...,\widehat{Y}_{n})=\widehat{\kappa}_{n}\left(\left(\frac{X_{1}}{\widehat{c}_{n}}\right)^{\widehat{\kappa}_{n}},...,\left(\frac{X_{n}}{\widehat{c}_{n}}\right)^{\widehat{\kappa}_{n}}\right)=1.$

###### Remark 3.2

The invariance properties figuring in (3.1) and (3.2) are not restricted to
the MLE alone. Other estimators, such as moment estimators or estimators based
on order statistics, may also result in test statistics that satisfy these
properties provided that the estimators under discussion satisfy the
equivariance properties referred to in the Introduction. A particular case of
equivariant estimators with minimum risk are the Pitman estimators of location
and scale; see Zacks (1971)[§7.2], Lehmann and Casella (1998)[§3.1] and
Lin’kov (2005)[§3.1-3.2].

In this connection, and before closing this section we wish to emphasize that
the invariance properties figuring in (3.1), (3.2), and (3.6), do not imply
that our test procedures reduce to the case of simple hypotheses with
corresponding parameters known. In fact parameter estimation generally does
have an effect on the distributional properties of the tests, and test
invariance only means that these distributional properties do not involve the
actual true values of the unknown parameters being estimated; see for instance
Meintanis and Swanepoel (2007). In the next section we illustrate the
performance of some GOF tests for exponentiality that are employed in order to
test for the Weibull distribution with both parameters unknown.

##  4 Monte Carlo and real–data

In this section we study the finite–sample performance of a few tests for the
Weibull distribution with DF, $F_{c,\kappa}(x)=1-\exp\\{-(x/c)^{\kappa}\\}$,
and unknown parameter $(c,\kappa)$. Recall that if the tests are applied on
$(\widehat{Y}_{j},\ j=1,\ldots,n)$ as defined in (3.5) with
$(\widehat{c}_{n},\widehat{\kappa}_{n})$ being the MLE, then we can set
$c=\kappa=1$, and consequently we may invoke any test for unit exponentiality
(see item 1 in Section 1). From the plethora of available tests we consider
the Anderson–Darling test based on

${\rm{AD}}_{n}=-n-\frac{1}{n}\sum_{j=1}^{n}(2j-1)(\log{Z_{(j)}}+\log(1-Z_{(n+1-j)})),$

where $Z_{(j)}=1-{\rm{e}}^{-\widehat{Y}_{(j)}}\ (j=1,...,n)$ with
$\widehat{Y}_{(1)}\leq\widehat{Y}_{(2)}{\leq}\ldots\leq\widehat{Y}_{(n)}$
being the ordered statistics. The ADn test is often the most powerful test
among the classical tests based on the empirical DF. We also include the test
of Henze and Meintanis (2002) based on

${\rm{HM}}_{n}=\frac{1}{n}\sum_{j,k=1}^{n}\frac{1+(\widehat{Y}_{j}+\widehat{Y}_{k}+2)^{2}}{(\widehat{Y}_{j}+\widehat{Y}_{k}+1)^{3}}-2\sum_{j=1}^{n}\frac{\widehat{Y}_{j}+2}{(\widehat{Y}_{j}+1)^{2}}+n,$

which is amongst the best performing exponentiality tests in the comparison
studies of Henze and Meintanis (2005) and Allison et al. (2017), and a smooth
test of fit for the exponential distribution (see Rayner et al. (2009)[§6.3]
given by

${\rm{RB}}_{n}=\frac{1}{n}\left(\sum_{j=1}^{n}L_{2}(\widehat{Y}_{j})\right)^{2}+\frac{1}{n}\left(\sum_{j=1}^{n}L_{3}(\widehat{Y}_{j})\right)^{2},$

where

$L_{2}(z)=1-2z+z^{2}/2,\ L_{3}(z)=1-3z+3z^{2}/2-z^{3}/6,$

are the Laguerre polynomials of orders 2 and 3.

We consider tests of size $\alpha$, and for a given sample size $n$ we
calculate the test statistics based on a large number $M$ of Monte Carlo
samples and obtain the quantile corresponding to $1-\alpha$. Specifically for
each sample of size $n$ and each test statistic, say $T$, we generate
observations from a Weibull distribution with a fixed combination of
$(c,\kappa)$, then we calculate the MLE
$(\widehat{c}_{n},\widehat{\kappa}_{n})$ and obtain the value of the test
statistic $T_{m}$ based on the transformed sample
$\widehat{Y}_{j}=(X_{j}/\widehat{c}_{n})^{\widehat{\kappa}_{n}},\ j=1,...,n$.
By iterating this procedure for $m=1,2,...,M$, we obtain the critical value of
the test statistic as the $1-\alpha$ quantile of the empirical distribution of
$(T_{m},\ m=1,...,M)$. By the test invariance articulated in the previous
section we only need to draw samples from the Weibull distribution with
$(c,\kappa)=(1,1)$, i.e. by sampling from the unit exponential distribution.
Nevertheless we examined the three test statistics by sampling from a Weibull
distribution with several combinations of $(c,\kappa)$ and indeed our
conclusion for a parameter–free test statistic was confirmed as the resulting
critical values remained stable regardless of the actual value of $(c,\kappa)$
employed.

The actual Monte Carlo was performed with sample size $n=50,100,150$ and
$n=200$, with $M=100,000$ iterations at significance level $\alpha=0.1,0.05$
and $\alpha=0.01$, and the resulting critical values are reported in Table
LABEL:table:criticalpt. The figures in Table LABEL:table:criticalpt show that
convergence to the asymptotic distribution is faster for the ADn and HMn
tests, while the smooth test RBn is somewhat slower to reach its limit
distribution.

Table 1: Critical values for ADn, HMn, and RBn, at significance level $\alpha$ | ADn | HMn | RBn
---|---|---|---
$n\setminus\alpha$ | $0.1$ | $0.05$ | $0.01$ | $0.1$ | $0.05$ | $0.01$ | $0.1$ | $0.05$ | $0.01$
$n=50$ | $0.629$ | $0.750$ | $1.027$ | $0.036$ | $0.047$ | $0.076$ | $0.922$ | $1.374$ | $4.558$
$n=100$ | $0.623$ | $0.755$ | $1.027$ | $0.037$ | $0.048$ | $0.078$ | $1.123$ | $1.797$ | $7.456$
$n=150$ | $0.634$ | $0.755$ | $1.041$ | $0.037$ | $0.049$ | $0.079$ | $1.269$ | $2.092$ | $9.249$
$n=200$ | $0.634$ | $0.757$ | $1.050$ | $0.038$ | $0.049$ | $0.078$ | $1.383$ | $2.313$ | $9.914$

Using this methodology we apply the three tests on two real–data sets of sizes
$n=63$ (Ex.1) and $n=46$ (Ex.2) employed by Smith and Naylor (1987) (see Table
1, p.359 for the data). The data correspond to experimental measurements on
the strength of glass fiber of length $1.5$ cm (Ex.1) and $15$ cm (Ex.2).

Table 2: Actual values of test statistics (“act.”), and critical values at significance level 0.1, 0.05, 0.01 | ADn | HMn | RBn
---|---|---|---
| act. | $0.1$ | $0.05$ | $0.01$ | act. | $0.1$ | $0.05$ | $0.01$ | act. | $0.1$ | $0.05$ | $0.01$
Ex.1 | $\underline{1.241}$ | $0.629$ | $0.747$ | $1.031$ | $\underline{0.104}$ | $0.036$ | $0.047$ | $0.076$ | $\underline{1.458}$ | $0.981$ | $1.496$ | $5.820$
Ex.2 | $\underline{0.329}$ | $0.631$ | $0.752$ | $1.024$ | $\underline{0.021}$ | $0.036$ | $0.046$ | $0.077$ | $\underline{0.602}$ | $0.904$ | $1.349$ | $4.277$

The values of the MLE are
$(\widehat{c}_{n},\widehat{\kappa}_{n})=(5.781,1.628)$ for Ex.1, and
$(\widehat{c}_{n},\widehat{\kappa}_{n})=(5.147,1.230)$ for Ex.2. From these
estimated values, we obtain $(\widehat{Y}_{j},\ j=1,...,n)$, and calculate
each of the three test statistics, whose actual values are reported in Table
LABEL:table:actual (underlined figures). Before further analysis we also
considered the Kolmogorov-Smirnov (KS) test with the data $\widehat{Y}_{j}$,
and thereby obtained p-values $0.1078$ for Ex.1 and $0.9473$ for Ex.2. These
are consistent with the values reported in Table 4 (Ex.1) and Table 5 (Ex.2),
of the real–data analysis in §5 of Wu et al. (2021). Thus at significance
level $10\%$, and on the basis of the KS test we can not reject the hypothesis
of exponentiality of the standardized data $(\widehat{Y}_{j},\ j=1,...,n)$,
which in turn can be interpreted to imply that the source data $(X_{j},\
j=1,...,n)$ might have originated from a Weibull distribution with the MLE
estimates as parameters. On the other hand, the exponentiality of the data of
Ex.1 is rejected by the ADn and the HMn tests, and only the smooth test at
significance levels $5\%$ and $1\%$ finds no evidence to reject this
hypothesis. At the same time, the exponentiality for the data of Ex.2 is
supported by all three tests uniformly over the values of $\alpha$ considered.
The corresponding histogram and distribution function plots that are shown in
Figure 1 further corroborate our results. Therefore, there is strong evidence
in favour of a Weibull distribution for the data of Ex.2, while the
corresponding conclusion for the data of Ex.1 should be questioned and
occasional non–rejection may be due to low power. For instance, the KS test is
often the least powerful amongst the classical tests based on the empirical
DF.

Figure 1: Histogram and distribution function of the standardized strength
data of glass fiber (Ex.1 and Ex.2). The standard exponential density and
distribution function are superimposed on the corresponding graphs (solid
curved lines).

## References

  * Allison et al. (2017) Allison, J.S., Santana, L., Smit, N. and Visagie, I.J.H. (2017) An ‘apples to apples’ comparison of various tests for exponentiality. Computational Statistics 32, 1241–1283.
  * Drost et al. (1990) Drost, F.C., Kallenberg, W.C.M. and Oosterhoff, J. (1990) The power of EDF tests of fit under nonrobust estimation of nuisance parameters. Statistics & Decisions 8, 167–182.
  * Epps (1993) Epps, T.W. (1993) Characteristic functions and their empirical counterparts: Geometrical interpretations and applications to statistical inference. The American Statistician 47, 33–38.
  * Epps (2005) Epps, T.W. (2005) Tests for location-scale families based on the empirical characteristic function. Metrika 62, 99–114.
  * Epps and Pulley (1983) Epps, T.W. and Pulley, L.B. (1983) A test for normality based on the empirical characteristic function. Biometrika 70, 723–726.
  * Escanciano (2009) Escanciano, J.C. (2009) On the lack of power of omnibus specification tests. Econometric Theory 25, 162–194.
  * Hall et al. (2013) Hall, P., Lombard, F. and Potgieter, C.J. (2013) A new approach to function-based hypothesis testing in location-scale families. Technometrics 55, 215–223.
  * Hallin and Jurečková (2012) Hallin, M. and Jurečková, J. (2012) Equivariant estimation. In: Encyclopedia of Environmetrics, 2nd Edition, DOI: 10.1002/9780470057339.vnn082. J. Wiley & Sons, New York.
  * Henze (2002) Henze, N. (2002) Invariant tests for multivariate normality: a critical review. Statistical Papers 43, 467–506.
  * Henze and Meintanis (2002) Henze, N. and Meintanis, S.G. (2002) Tests of fit for exponentiality based on the empirical Laplace transform. Statistics 36, 147–161.
  * Henze and Meintanis (2005) Henze, N. and Meintanis, S.G. (2005) Recent and classical tests for exponentiality: a partial review with comparisons. Metrika 61, 29–45.
  * Janssen (2000) Janssen, A. (2000) Global power functions of goodness of fit tests. The Annals of Statistics 28, 239–253.
  * Klar (2000) Klar, B. (2000) Diagnostic smooth tests of fit. Metrika 52, 237–252.
  * Ledwina (1994) Ledwina, T. (1994) Data–driven version of Neyman’s smooth test of fit. Journal of the American Statistical Association 89, 1000–1005.
  * Lehmann and Casella (1998) Lehmann, E.L. and Casella, G. (1998) Theory of Point Estimation. Springer, New York.
  * Lehmann and Romano (2005) Lehmann, E.L. and Romano, J.P. (2005) Testing Statistical Hypotheses. Springer, New York.
  * Lin’kov (2005) Lin’kov, Yu.N. (2005) Lectures in Mathematical Statistics Parts 1 and 2. American Mathematical Society, Rhode Island.
  * Meintanis et al. (2015) Meintanis, S.G., Ngatchou–Wandji, J. and Taufer, E. (2015) Goodness-of-fit tests for multivariate stable distributions based on the empirical characteristic function. Journal of Multivariate Analysis 140, 171–192.
  * Meintanis and Swanepoel (2007) Meintanis, S.G. and Swanepoel, J. (2007) Bootstrap goodness-of-fit tests with estimated parameters based on empirical transforms. Statistics & Probability Letters 77, 1004–1013.
  * Neyman (1937) Neyman, J. (1937) Smooth test for goodness of fit. Skandinavisk Aktuarietidskrift 20, 149–199.
  * Rayner et al. (2009) Rayner, J.C.W., Thas, O. and Best, D.J. (2009) Smooth Tests of Goodness of Fit: using R. Wiley & Sons, New York.
  * Sinner et al. (2023) Sinner, C., Dominicy, Y., Trufin, J., Waterschoot, W., Weber, P. and Ley, C. (2023) From Pareto to Weibull– A constructive review of distributions on $\mathbb{R}^{+}$. International Statistical Review 91, 35–54.
  * Smith and Naylor (1987) Smith, R.L. and Naylor, J.C. (1987) A comparison of maximum likelihood and Bayesian estimators for the three- parameter Weibull distribution. Journal of the Royal Statistical Society. Series C (Applied Statistics) 36, 358–369.
  * Thode (2002) Thode, H.C. (2002) Testing for Normality. Marcel Dekker, New York.
  * Vexler and Hutson (2023) Vexler, A. and Hutson, A.D. (2023) A characterization of most(more) powerful test statistics with simple nonparametric applications. The American Statistician, DOI: 10.1080/00031305.2023.2192746.
  * Wu et al. (2021) Wu, Y., Xie, H., Chiang, J.Y., Peng, G. and Qin, Y. (2021) Parameter estimation and applications of the Weibull distribution for strength data of glass fiber. Mathematical Problems in Engineering, DOI: 10.1155/2021/9175170.
  * Zacks (1971) Zacks, S. (1971) The Theory of Statistical Inference. John Wiley & Sons, New York.


# Machine-learning applied to classify flow-induced sound parameters from
simulated human voice

Florian Kraxberger
Institute of Fundamentals and Theory in Electrical Engineering (IGTE)
Graz University of Technology
8010 Graz, Austria
<EMAIL_ADDRESS>
& Andreas Wurzinger
Institute of Fundamentals and Theory in Electrical Engineering (IGTE)
Graz University of Technology
8010 Graz, Austria
<EMAIL_ADDRESS>
& Stefan Schoder
Institute of Fundamentals and Theory in Electrical Engineering (IGTE)
Graz University of Technology
8010 Graz, Austria
<EMAIL_ADDRESS>

###### Abstract

Disorders of voice production have severe effects on the quality of life of
the affected individuals. A simulation approach is used to investigate the
cause-effect chain in voice production showing typical characteristics of
voice such as sub-glottal pressure and of functional voice disorders as
glottal closure insufficiency and left-right asymmetry. Therewith, 24
different voice configurations are simulated in a parameter study using a
previously published hybrid aeroacoustic simulation model. Based on these 24
simulation configurations, selected acoustic parameters (HNR, CPP, …) at
simulation evaluation points are correlated with these simulation
configuration details to derive characteristic insight in the flow-induced
sound generation of human phonation based on simulation results. Recently,
several institutions studied experimental data, of flow and acoustic
properties and correlated it with healthy and disordered voice signals. Upon
this, the study is a next step towards a detailed dataset definition, the
dataset is small, but the definition of relevant characteristics are precise
based on the existing simulation methodology of simVoice. The small datasets
are studied by correlation analysis, and a Support Vector Machine classifier
with RBF kernel is used to classify the representations. With the use of
Linear Discriminant Analysis the dimensions of the individual studies are
visualized. This allows to draw correlations and determine the most important
features evaluated from the acoustic signals in front of the mouth. The GC
type can be best discriminated based on CPP and boxplot visualizations.
Furthermore and using the LDA-dimensionality-reduced feature space, one can
best classify subglottal pressure with 91.7% accuracy, independent of healthy
or disordered voice simulation parameters.

_Keywords_ Voice Signal $\cdot$ Voice Disorders $\cdot$ Aeroacoustics $\cdot$
Phonation

##  1 Introduction

Human phonation is the result of complex aeroacoustic phenomena such as
periodical interruptions of the glottal airstream resulting in a pulsatile jet
flow in the larynx, which are the main source of the human voice [1]. Using
numerical models enable to analyze the cause-effect chains in the complex
human phonation system. Thereby, for the driven vocal folds simulation model
[2] parameters constitute the ’cause’ (i.e. regular and irregular oscillation
configurations, sub-glottal pressure, …), and the simulation results the
’effect’ (i.e. flow field, acoustic pressure signals, … ). According to [1],
the efficiency of human phonation is mainly influenced by the vocal fold
closure characteristics and the vocal fold motion properties. Thus, a high
phonation efficiency is present, when the vocal folds close the glottis
completely in each oscillation cycle, and when the vocal folds oscillate
symmetrically as well as periodically. The phonation quality is reduced, when
the glottis does not close completely or oscillates asymmetrically, resulting
in a breathy or hoarse voice quality due to increased broadband components in
the voice signal [3, 4, 5]. These irregularities can occur in visually healthy
larynges [6, 7, 8] as well as pathologic larynges [9] and their occurrence
correlates with increasing patient age [10, 11].

Recently machine learning-based approaches improved the insights in voice
research [12, 13, 14] and to separate healthy from disordered voices [15, 16,
17]. Using a self-organizing map based on acoustic parameters, [15]
differentiated normal from disordered voices with 0.7628 accuracy. In [16] and
for separation of normal, breathy, rough and hoarse voices using acoustic
parameters 0.75 accuracy was achieved. Higher accuracy of 0.81, was achieved
in [17] using PVG based parameters. Also, other studies showed good prediction
capabilities by using acoustic measures [18, 19]. Somewhat the pure
differentiation between healthy and some features associated with functional
dysphonia (for instance incomplete glottis closure) is not precise since they
also occur in healthy subjects [20] and some parameters are redundant [21,
22]. Regarding the broad range of possible impacts on the voice, experimental
study data were enhanced with machine learning techniques to get a more
profound understanding about the redundancy of which parameters best
characterize functional dysphonia or other voice disorders [23].

Differently to experimental approaches (also using machine learning) and since
disorders of voice production have severe effects on the quality of life, a
simulation model simVoice [24] was developed to investigate the cause-effect
chain in voice production with fully resolved data. The simVoice model depicts
typical characteristics of voice such as sub-glottal pressure and of
functional voice disorders as glottal closure insufficiency and left-right
asymmetry. In total up to 80 different voice configurations were simulated.
This type of simulation results in high-resolution data that can be used for
machine learning techniques as input. The motivation of this paper is to
demonstrate the usability of and to perform systematic investigation on the
voice signal (similar to [19]). The following simulation input parameters were
used: (i) four different levels of glottis closure, (ii) symmetric and
asymmetric vocal fold motion of the same frequency, and (iii) three different
subglottal pressures, to analyze the relation between voice signal measures
such as HNR, CPP and more. Also, we have to note that this is a first study on
using simulation results for machine learning in the human voice community.
Furthermore, the result interpretation is limited due to the small dataset.
Nevertheless, interesting conclusions can be drawn from this high-quality,
high-resolution simulation dataset.

##  2 Methods

In previous works, a validated hybrid aeroacoustic model of human phonation
has been established [2, 25, 26, 27], which is called simVoice. To reduce the
computational cost of a fluid-structure-acoustic interaction model [27], a
prescribed vocal fold motion model is used [28]. This model consists of a
finite-volume incompressible fluid dynamics simulation, and a finite-element
acoustic simulation, as documented in [29]. In [24], the model has been used
to model different regular and irregular phonation characteristics, such as
various degrees of insufficient glottal closure, symmetric or asymmetric vocal
fold motion, and different subglottal pressures. This section is a brief
introduction to the numerical model simVoice and the acoustic features of the
simulated voice signals are defined. The aim is to discriminate against
different healthy and functional-based voice disorders solely based on the
evaluated acoustic features (e.g. [19]). This approach originates in voice and
music research and usually uses a large number (often
$N_{\text{features}}>20$) of audio features [30, 31, 32, 33]. Regarding the
small dataset and the often redundant features [21], only a few features were
selected for the present study are: (i) sound pressure level (SPL), (ii)
Harmonics-to-Noise Ratio (HNR), (iii) Cepstral Peak Prominence (CPP), (iv)
spectral slope, (v) Hammarberg Index (HBI), and (vi) Alpha Ratio (see Sec.
2.3).

###  2.1 Hybrid Aeroacoustic Simulation Setup for Human Phonation (simVoice)

The hybrid aeroacoustic simulation setup used to model human phonation
consists of four components: (i) an incompressible computational fluid
dynamics (CFD) simulation, (ii) conservative interpolation from the dense grid
of the CFD to the coarser grid of the computational acoustics (CA) simulation,
(iii) computation of the aeroacoustic source term and (iv) acoustics
simulation based on the Perturbed Convective Wave Equation (PCWE).

The incompressible CFD simulation is documented in [26, 28, 27] and has been
performed with the software StarCCM+. The result thereof is the incompressible
pressure $p^{\mathrm{ic}}$, as described in [24], and [25] provides an in-
depth analysis of the aeroacoustic source terms. The aeroacoustic source terms
are furthermore filtered with the algorithm proposed in [34]. Figure 1 shows
the computational domain of the CFD simulation.

Figure 1: (color online) Computational domain of the CFD. The subglottal
pressure is applied as a constant inlet pressure boundary condition [24, Fig.
2(A)].

The PCWE, as defined in [35, 2], is a wave equation regarding the acoustic
scalar potential $\psi^{\mathrm{a}}$. From the acoustic scalar potential
$\psi^{\mathrm{a}}$, the acoustic pressure $p^{\mathrm{a}}$ is evaluated.

It has been shown in [25], that the convective term can be neglected in the
substantial derivative, without a significant deviation in the spectrum of the
resulting acoustic pressure. The PCWE is solved using the finite-element
solver openCFS [36]. An in-depth description of the CA simulations is
available in [29, 2]. Fig. 2 shows the domain of the CA simulations consisting
of the larynx, vocal tract, propagation and perfectly matched layer (PML)
regions. The CA simulation domain presents an extension of the CFD simulation
domain shown in Fig. 1 with a propagation and PML region. For the inlet
surface, an absorbing boundary condition (ABC) was used in the CA simulations.
The vocal tract and propagation regions are connected with a non-conforming
(NC) interface, as described in [29]. The PML region ensures free-field wave
propagation for the propagation domain [37]. The walls of the vocal tract are
modeled sound-hard.

(a) Slice of the CA domain

(b) CA domain

Figure 2: (color online) Computational domain of the CA simulations.

###  2.2 Phonation configurations (labels)

Voice disorders are phenomena containing complex interactions between
different parts of the voice production process, such as the subglottal
pressure coming from the bronchial tubes, the vocal folds (oscillating at a
frequency of $148\text{\,}\mathrm{Hz}$ representing male voice), the false
vocal folds, the vocal tract (/a/ [38]). For the present study, these
interacting elements and regular (subglottal pressure) and irregular
characteristics (asymmetry and insufficiency) thereof are abstracted in terms
of certain functions of the voice production process (see [24]). Thus,
periodic symmetric and asymmetric oscillation characteristics, different
grades of glottis insufficiency, in combination with varying subglottal
pressure have been simulated as described in [24, 39]. The simulation results
of the acoustic pressure signal in front of the vocal tract (the mouth), are
used for the later computations of the acoustic features.

##### Subglottal Pressure

The presence of a pressure difference between the subglottal and supraglottal
regions drives the glottal airstream. The inability to build up such a
pressure can be one reason for a dysphonic voice, e.g. caused by muscle
tension disorders, as reported in [40, 41]. The simulation cases include the
subglottal pressures as follows

$P_{\text{inlet}}\in\left\\{385,775,1500\right\\}$\text{\,}\mathrm{Pa}$.$ (1)

Where, $775\text{\,}\mathrm{Pa}$ represents normal subglottal pressure and the
other low or high, respectively.

##### Glottal Closure Type

The vocal folds can have varying degrees of glottal insufficiency, which occur
when the vocal folds do not close completely in each cycle. According to [42,
43], this irregularity occurs for many voice pathologies. As a consequence,
concerned patients report, that a higher effort is necessary for phonation
[44]. The glottal closure (GC) type determines the amount of glottal
insufficiency. $o_{\mathrm{initial}}$ is defined as the percentage of the
vocal folds lengths which does not close. The simulation model employs initial
glottal openings of

$o_{\mathrm{initial}}\in\left\\{0\text{(GC1)},40,70,100\text{(GC4)}\right\\}$\text{\,}\mathrm{\char
37\relax}$,$ (2)

which are called GC1 (complete glottal closure), GC2 (60% of the vocal folds
lengths close), GC3 (30% of the vocal folds lengths close), and GC4 (no
contact of vocal folds during oscillations), respectively [24]. The different
GC types are based on high-speed videos from experiments with excised larynges
and are motivated by clinical observations [45, 46, 27]. A detailed
description as well as visualizations of the different symmetry cases are
available in Fig. 5.5 [47].

##### Vocal Fold Motion Symmetry

In regular phonation, the left and right vocal fold oscillate symmetrically.
However, in the case of laryngeal unilateral paresis or unbalanced muscle
tension, the vocal folds oscillate asymmetrically [48, 49]. Therefore, the
proposed simulation model features two cases of vocal fold motion symmetry. In
the symmetric case, both vocal folds have the same amplitude in their motion.
In the asymmetric case, one vocal fold oscillates with
$50\text{\,}\mathrm{\char 37\relax}$ of the amplitude of the other vocal fold.
Both symmetry cases are depicted in Fig. 3.

Figure 3: (color online) Exemplary vocal fold motion of GC1 for the symmetric
and asymmetric case along the y-axis (medial-lateral direction) for the center
point on the medial plane of the VF surface [24, Fig. 4].

In Tab. 1, the simulation parameters are summarized. A total of 24 different
functional-based voice disorders have been simulated.

Table 1: Simulation parameters of functional-based voice disorders as defined in [24]. Parameter | # | Values
---|---|---
Subglottal Pressure | 3 | $P_{\text{inlet}}\in\left\\{385,775,1500\right\\}$\text{\,}\mathrm{Pa}$$
Glottal Closure Type | 4 | $o_{\mathrm{initial}}\in\left\\{0,40,70,100\right\\}$\text{\,}\mathrm{\char 37\relax}$$
Vocal Fold Motion Symmetry | 2 | asymmetric or symmetric
| 24 | configurations total

###  2.3 Acoustic Feature Definitions

The acoustic pressure signal results of the simVoice model are recorded at an
evaluation point (so-called ’microphone’), as depicted in Fig. 2 (a). The
recorded signal is time discrete acoustic pressure $p[n]$ with $n$ being the
time index [50] sampled at
$T_{s}=1/f_{s}=$13.6\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ [2, 24]. With
those audio signals (one resolved signal for each of the 24 simulation
configurations), we compute certain acoustic features for the analysis. Thus,
we aim to analyze the flow-induced sound generation of health and of
functional-based voice disorders by using the evaluated features only.

##### Sound Pressure Level

The SPL $L_{p}$ of an acoustic pressure signal $p[n]$ is computed as follows,

$\begin{split}L_{p}&=20\cdot\log_{10}{\frac{\tilde{p}}{p^{\mathrm{a}}_{\mathrm{ref}}}}\qquad\text{with}\qquad
p^{\mathrm{a}}_{\mathrm{ref}}=$20\text{\,}\mathrm{\SIUnitSymbolMicro Pa}$\\\
\tilde{p}&=\sqrt{\frac{1}{N}\sum_{n=1}^{N}\left(p[n]\right)^{2}}.\end{split}$
(3)

Thereby, $N$ is the number of time steps of $p[n]$. Consequently, $\tilde{p}$
denotes the root-mean-square of the acoustic pressure signal. The subglottal
pressure is expected to have a significant influence on the sound pressure
level (SPL) magnitude of the acoustic signal. Since a higher subglottal
pressure (driving force of the phonation) will certainly exhibit a stronger
voice at the same vocal efficiency.

##### Harmonics-to-Noise Ratio

Harmonics-to-Noise ratio (HNR) provides a measure for the degree of hoarseness
[51], and has since been used, among others, to identify voice disorders [52]
and voice quality [53, 54]. Thus, in the present study HNR is used with the
aim of contributing to the identification and analysis of glottis closure
type, the asymmetry and their combinations defined in Sec. 2.2. The idea is to
use the HNR to be a measure of an increasing amount of turbulent structures
and decreasing tonal contribution during the whole oscillation cycles for only
partly closing vocal folds. HNR is defined as the energy ratio between the
harmonic signal component to the noise-like signal component in
$\text{\,}\mathrm{dB}$. The splitting of harmonic and noise-like components is
performed via the auto-correlation function (ACF), as described in [31, pp.
77–78]. The ACF $r_{p}(\tau)$ of a pressure signal $p[n]$ is calculated with

$r_{p}(\tau)=\sum_{n=0}^{N-1}p[n]p[n+\tau],$ (4)

where $\tau$ is the independent lag variable. Therewith, the HNR can be
computed with [31, Eq. (2.234)], such that

$\mathrm{HNR}_{\mathrm{dB,Eyben}}=10\log_{10}{\frac{\max\left\\{r_{p}(\tau)\right\\}}{r_{p}(0)-\max\left\\{r_{p}(\tau)\right\\}}},$
(5)

under the assumption that (i) the noise is additive, (ii) it is uncorrelated
with the signal and (iii) the noise is uncorrelated to itself (i.e. it is
white noise). To be robust against deviations from the assumptions in the
simulated signal, a margin of $k_{\mathrm{marg}}=300$ lag bins is introduced,
which corresponds to a lag of
$\tau_{\mathrm{marg}}=k_{\mathrm{marg}}/f_{s}=$4.08\text{\,}\mathrm{ms}$$.
Therewith, the HNR $h$ is calculated as follows:

$\mathrm{HNR}_{\mathrm{dB}}=h=10\log_{10}{\frac{\max_{\tau_{\mathrm{marg}}}^{\tau_{\mathrm{max}}}\left\\{r_{p}(\tau)\right\\}}{r_{p}(0)-\max_{\tau_{\mathrm{marg}}}^{\tau_{\mathrm{max}}}\left\\{r_{p}(\tau)\right\\}}},$
(6)

where $\tau_{\mathrm{max}}$ is the maximum possible lag time (i.e. the time
corresponding to the highest lag index).

##### Cepstral Peak Prominence

The harmonic structures occurring in the voice spectrum map to peaks in the
cepstral domain, as introduced in [55]. Healthy voice shows, that the harmonic
parts dominate the spectrum, and thus it is expected that the corresponding
peak in the cepstral domain is prominent compared to non-harmonic parts of the
cepstrum. Hence, Cepstral Peak Prominence (CPP) measures the prominence of the
peak in the cepstral domain and therefore the prominence of harmonic
structures in the spectrum compared to non-harmonic parts. Therewith, we aim
to discriminate simulation parameters which affect the harmonic structure of
the voice signal, e.g. by an incomplete glottis closure. The CPP is computed
as follows [56, 57]: Given a pressure signal $p[n]$, a power spectral density
(PSD) estimate $\tilde{p}[k]$ is computed using Welch’s method [58] such that

$\begin{split}\tilde{p}[k]&=\sqrt{\mathcal{F}_{n\mapsto
k}\\!\left\\{p[n]\right\\}\\!\left[k\right]}\\\
\tilde{p}_{\mathrm{dB}}[k]&=20\cdot\log_{10}{\frac{\tilde{p}[k]}{p_{0}}}\qquad\text{with}\qquad
p_{0}=$20\text{\,}\mathrm{\SIUnitSymbolMicro Pa}$,\end{split}$ (7)

where $\mathcal{F}_{n\mapsto k}\\!\left\\{p[n]\right\\}\\!\left[k\right]$
denotes the PSD based on the $N$-point discrete Fourier transform (DFT) [50].
In accordance to [55], the cepstrum $c_{\tilde{p}}[k]$ is computed such that

$c_{\tilde{p},\mathrm{dB}}[q]=10\cdot\log_{10}{\left|\mathcal{F}_{k\mapsto
q}\\!\left\\{\tilde{p}_{\mathrm{dB}}[k]\right\\}\\!\left[q\right]\right|^{2}}.$
(8)

The independent variable of the cepstrum is called quefrency and has the
dimension of time [55]. The CPP $c$ is evaluated as follows [56, 57]: For
quefrencies larger than $1\text{\,}\mathrm{ms}$, a linear regression line is
fitted to the cepstrum. The CPP is the difference between the peak amplitude
of the cepstral maximum and the linear regression line evaluated at the
quefrency of the cepstral maximum. Thus, it represents how far the cepstral
peak emerges from the cepstral "background noise". This can be interpreted as
the prominence of harmonics in the signal compared to non-harmonic signal
portions. In [46, 59], it was used together with other features as a measure
for vocal effort and to disambiguate different voice disorders. In [60, 61,
62], CPP and HNR are used to compute the so-called Acoustic Voice Quality
Index together with shimmer and spectral slope and tilt. Furthermore, CPP and
cepstrum-based features are used as fundamental measures in speech sciences
[63, 64, 65].

##### Spectral Slope

The spectral slope is a measure of energy distribution across the spectrum and
is used in [66, 67] to distinguish voice qualities. If the glottal airstream
is interrupted off completely in each vocal fold oscillation cycle, we expect
effects on the energy distribution in the resulting acoustic signal. Thus,
this feature is used to determine if high frequencies (which we expect to
increase due to incomplete glottis closures) dominate lower frequencies, such
that the general spectral trend represented by the spectral slope is
influenced. In [31], spectral slope111As shown in [68], the terms ”spectral
slope” and ”spectral tilt” are used interchangeably in the literature, thus,
for spectral slope, the definition of [31, pp. 35] is used. is defined as the
result of a linear regression operation on the whole spectrum $\tilde{p}[k]$.
Thereby, a line $\hat{y}=ax+b$ is fitted to the spectrum
$\tilde{p}[k]=\sqrt{\mathcal{F}_{n\mapsto
k}\\!\left\\{p[n]\right\\}\\!\left[k\right]}$, such that the mean-squared
error between the line and the original spectrum is minimized (minimum mean-
squared error condition). The parameter of interest is $a$, which is the slope
of the linear regression result. A formula for $a$ is provided by [31], such
that

$a=\frac{N\Sigma_{xy}-\Sigma_{x}\Sigma_{y}}{N\Sigma_{x^{2}}-\Sigma_{x}^{2}},$
(9)

where $N$ is the number of DFT frequency bins and $\Sigma_{xy}$, $\Sigma_{x}$,
$\Sigma_{y}$ and $\Sigma_{x^{2}}$ are defined as follows

$\begin{split}\Sigma_{x}=\sum_{i=0}^{N-1}f(k_{i}),&\qquad\Sigma_{y}=\sum_{i=0}^{N-1}p[i],\\\
\Sigma_{x^{2}}=\sum_{i=0}^{N-1}f(k_{i})^{2},&\qquad\Sigma_{xy}=\sum_{i=0}^{N-1}f(k_{i})p[i],\end{split}$
(10)

where $f(k)$ is the frequency in $\text{\,}\mathrm{Hz}$ associated with the
frequency bin $k$. By exploiting the equidistant frequency scale used in the
DFT, simplifications arise [31], such that

$\Sigma_{x}=\frac{1}{2}N(N-1)\qquad\text{and}\qquad\Sigma_{x^{2}}=\frac{1}{6}N(N-1)(2N-1).$
(11)

##### Hammarberg Index and Alpha Ratio

Hammarberg Index (HBI) and Alpha Ratio are measures of the energy distribution
across the spectrum. Similar to spectral slope, also these features are
introduced aiming to discriminate simulation parameters which influence the
glottis closure, i.e. GC type and glottis symmetry. HBI $\eta$ as defined in
[69] is the ratio of the most prominent energy peak in the range of
$0\text{\,}\mathrm{kHz}$–$2\text{\,}\mathrm{kHz}$ to the most prominent energy
peak in the region $2\text{\,}\mathrm{kHz}$–$5\text{\,}\mathrm{kHz}$. It is
calculated according to [31, p. 38], such that

$\eta=\frac{\max_{k=1}^{k_{\mathrm{pivot}}}\left\\{\tilde{p}[k]\right\\}}{\max_{k=k_{\mathrm{pivot}}+1}^{k_{\mathrm{max}}}\left\\{\tilde{p}[k]\right\\}},$
(12)

where the pivot frequency bin $k_{\mathrm{pivot}}$ is the highest spectral bin
for which $f\leq$2\text{\,}\mathrm{kHz}$$, and $k_{\mathrm{max}}$ is the
highest spectral bin for which $f\leq$5\text{\,}\mathrm{kHz}$$.

The alpha ratio is similar to the HBI described above, but instead of
computing the ratio between energy peaks, the energy sum in the frequency
bands is considered, as introduced in [70]. Using the frequency bands
$50\text{\,}\mathrm{Hz}$–$1\text{\,}\mathrm{kHz}$ and
$1\text{\,}\mathrm{kHz}$–$5\text{\,}\mathrm{kHz}$, the alpha ratio
$\rho_{\alpha}$ can be computed according to [31, pp. 38–39], such that

$\rho_{\alpha}=\frac{\sum\limits_{k=k_{\mathrm{start}}}^{k_{\mathrm{pivot}}}\tilde{p}[k]}{\sum\limits_{k=k_{\mathrm{pivot}}+1}^{k_{\mathrm{max}}}\tilde{p}[k]},$
(13)

where the $k_{\mathrm{start}}$ is the lowest spectral bin for which
$f\geq$50\text{\,}\mathrm{Hz}$$, $k_{\mathrm{pivot}}$ is the highest spectral
bin for which $f\leq$1\text{\,}\mathrm{kHz}$$, and $k_{\mathrm{max}}$ is the
highest spectral bin for which $f\leq$5\text{\,}\mathrm{kHz}$$.

##### Low-Pass Filtering

The acoustic pressure signals from the simVoice model is filtered by low-pass
filter. Before the computation of all features a filter with
$5\text{\,}\mathrm{kHz}$ cut-off frequency is used. This is due to the upper
frequency limit of the simulation setup [2, 34]. In addition, and before
computing CPP, HNR and spectral slope, we use a filter of
$2\text{\,}\mathrm{kHz}$ cut-off frequency targeting only the lowest two
formants [38]. The features based on both filtering strategies in total amount
to 9 evaluated acoustic features per configuration $i$. The evaluated features
for the $i$-th configuration ($i=1,\dots,24$) are denoted as a vector
$\bm{F}^{(i)}$, such that

$\bm{F}^{(i)}=\begin{bmatrix}L^{(i)}_{P,$5\text{\,}\mathrm{kHz}$}&h^{(i)}_{$5\text{\,}\mathrm{kHz}$}&h^{(i)}_{$2\text{\,}\mathrm{kHz}$}&c^{(i)}_{$5\text{\,}\mathrm{kHz}$}&c^{(i)}_{$2\text{\,}\mathrm{kHz}$}&a^{(i)}_{$5\text{\,}\mathrm{kHz}$}&a^{(i)}_{$2\text{\,}\mathrm{kHz}$}&\eta^{(i)}_{$5\text{\,}\mathrm{kHz}$}&\rho_{\alpha,$5\text{\,}\mathrm{kHz}$}^{(i)}\end{bmatrix}^{\mathrm{T}},$
(14)

where $L_{P}^{(i)}$ is SPL, $h^{(i)}$ is HNR, $c^{(i)}$ is CPP, $a^{(i)}$ is
spectral slope, $\eta^{(i)}$ is HBI, $\rho_{\alpha}^{(i)}$ is the alpha ratio,
and the subscript indicates the low-pass filter’s cut-off frequency.

###  2.4 Data Labeling, Dimensionality Reduction and Classification

We interpret different functional-based voice disorders as labels in a machine
learning sense. Thus, each of the 24 configurations defined in Sec. 2.2 is
labeled with one value of each of the following label sets. For the $i$-th
configuration, this yields labels from the following label sets:

$\begin{split}\ell^{(i)}_{\mathrm{subglottPress}}&\in\left\\{385,775,1500\right\\}\\\
\ell^{(i)}_{\mathrm{GC}}&\in\left\\{1,2,3,4\right\\}\\\
\ell^{(i)}_{\mathrm{sym}}&\in\left\\{0,1\right\\}.\end{split}$ (15)

Thus, the label $\bm{L}^{(i)}$ of the $i$-th simulation configuration is
formed as a vector

$\begin{split}\bm{L}^{(i)}=\begin{bmatrix}\ell^{(i)}_{\mathrm{subglottPress}}&\ell^{(i)}_{\mathrm{GC}}&\ell^{(i)}_{\mathrm{sym}}\end{bmatrix}^{\mathrm{T}}\end{split}$
(16)

For example, considering the configuration with $385\text{\,}\mathrm{Pa}$
subglottal pressure, GC1 and asymmetric VF motion, the label vector would be
$\bm{L}^{(1)}=\begin{bmatrix}385&1&0\end{bmatrix}^{\mathrm{T}}$. To form
classes, different oscillation patterns and subglottal pressure are considered
individually. Thereby, all evaluated feature vectors $\bm{F}^{(i)}$ are
considered to be of a certain class, if one element of the feature vector
$\bm{L}^{(i)}$ is identical. The other elements of the feature vector are
discarded. For example, the class GC1 contains all data samples
$\bm{F}^{(i)}$, for which
$\bm{L}^{(i)}=\begin{bmatrix}\ell^{(i)}_{\mathrm{subglottPress}}&1&\ell^{(i)}_{\mathrm{sym}}\end{bmatrix}^{\mathrm{T}}$,
regardless of the values of $\ell^{(i)}_{\mathrm{subglottPress}}$ and
$\ell^{(i)}_{\mathrm{sym}}$. The features, as defined in Sec. 2.3, span a
9-dimensional feature space. For dimension reduction and visualization, a two-
fold approach is used, consisting of (i) Linear Discriminant Analysis (LDA),
see Sec. 4.1.4 of [71] and (ii) upon the results of the LDA a classification
with Support Vector Machines (SVMs), see Sec. 7.1.3 of [71]. LDA is used to
project the data points $\bm{F}^{(i)}$ from the original feature space into a
two-dimensional feature space, with the goal to maximize the distance between
different classes while minimizing the covariance within each class. The
results of the LDA are transformed data points $\bm{F}^{(i)}_{\mathrm{tf}}$
for which $\mathrm{dim}\\!\left(\bm{F}^{(i)}_{\mathrm{tf}}\right)=2$ for the
LDA with respect to GC type (only using the first two major components for
visualization) and subglottal pressure, and
$\mathrm{dim}\\!\left(\bm{F}^{(i)}_{\mathrm{tf}}\right)=1$ for the LDA with
respect to symmetry property. In Fig. 4, the weights of the LDA transformation
are depicted.

Figure 4: Feature weights for dimension reduction.

Thereafter, the labeled data is used to train an SVM classifier with radial
basis function-kernels and a one-versus-rest classification. The SVM
implementation in the Python library sklearn has been used for this task [72].
To quantify the generalization ability of the SVM classifier, $5$-fold
stratified cross-validation was used [71]. The mean scores of the cross-
validation data sets are shown in Tab.LABEL:tab:03-acouPress-SVM-scores.
Thereafter, the whole data set was used to train the SVM classifier to
classify mesh grid points to visualize the class boundaries. Based on the
small number of data, the cross-validation procedure leads to stratified
results but should be considered in the future with increasing amount of data.

##  3 Results

First, we study Pearson’s correlation between the individual acoustic features
and display it as a correlation map. In Fig. 5, the correlation map is shown
for configurations evaluated from acoustic pressure at the evaluation position
as depicted in Fig. 2 (a).

Figure 5:  (color online) Correlation map between features and parameters for
signals (acoustic pressure at the microphone position). The fields are colored
according to the correlation value of the respective parameter combination
(i.e. the number in each field). The numbers in bold highlight cause-effect-
combinations which are further investigated using boxplot diagrams. The thick
black lines are separating simulation parameters from evaluated features.

From Fig. 5 it is visible, that for HNR, CPP and spectral slope, which have
been evaluated following a $2\text{\,}\mathrm{kHz}$ LP-filter provide stronger
correlations than those evaluated with a $5\text{\,}\mathrm{kHz}$ LP-filter.
Furthermore, it is visible that spectral slope does not correlate
significantly with any of the modeled phonation characteristics, i.e.
symmetry, GC type, and subglottal pressure. Based on the small dataset, the
significance of the individual correlation based on the autocorrelation is not
giving a full picture. Therefore, we continue studying the data with boxplots
to get a more detailed understanding of a significant difference of the data
points.

###  3.1 Boxplot Visualizations

In the following, boxplot diagrams are shown, which visualize the
distributions of the evaluated features. Thereby, manually selected 6 features
($L^{(i)}_{P,$5\text{\,}\mathrm{kHz}$},h^{(i)}_{$2\text{\,}\mathrm{kHz}$},c^{(i)}_{$2\text{\,}\mathrm{kHz}$},a^{(i)}_{$2\text{\,}\mathrm{kHz}$},\eta^{(i)}_{$5\text{\,}\mathrm{kHz}$},\rho_{\alpha,$5\text{\,}\mathrm{kHz}$}^{(i)}$)
are used because the correlation diagram in Fig. 5 already showed, which
features provide better correlations with respect to the simulation
parameters. From the boxplot diagrams, we can learn, how individual features
vary across different simulation parameters. They provide an insight on
possible parameter groups which cannot be discriminated using individual
features, but might be discriminable using a combination of features, e.g. by
means of an LDA/SVM approach.

##### Boxplots for Glottal Closure Type

In Fig. 6, the boxplot diagrams for the simulation parameter glottal closure
type are depicted for the acoustic pressure evaluated at the microphone
position. From Fig. 6 (A), (B), and (C), we see a general trend of decreasing
CPP and HNR values with increasing initial glottal openings. A significant
trend is visible for CPP evaluated with the low-pass filter at
$f_{c}=$2\text{\,}\mathrm{kHz}$$, as shown in Fig. 6 (B).

Figure 6: Boxplot diagrams of the features spectral slope $a$, CPP
$c_{$2\text{\,}\mathrm{kHz}$}$, HNR $h_{$2\text{\,}\mathrm{kHz}$}$, SPL
$L_{P,$5\text{\,}\mathrm{kHz}$}$, alpha ratio $\rho_{\alpha}$, and HBI $\eta$
as a function of glottal closure type for acoustic pressure signals at the
microphone position.

##### Boxplots for Symmetry

In Fig. 7, the boxplot diagrams for the simulation parameter symmetry are
depicted for the acoustic pressure evaluated at the microphone position. It is
shown, that no significant discrimination of glottal symmetry is possible
based on the evaluated features. This is also supported by near-zero
correlation values for the symmetry property depicted in Fig. 5. This non-
observable difference might be due to the inhomogeneous composition of the
data.

Figure 7: Boxplot diagrams of the features spectral slope
$a_{$2\text{\,}\mathrm{kHz}$}$, CPP $c_{$2\text{\,}\mathrm{kHz}$}$, HNR
$h_{$2\text{\,}\mathrm{kHz}$}$, SPL $L_{P,$5\text{\,}\mathrm{kHz}$}$, alpha
ratio $\rho_{\alpha}$, and HBI $\eta$ as a function of VF motion symmetry for
acoustic pressure signals at the microphone position.

##### Boxplots for Subglottal Pressure

In Fig. 8, the boxplot diagrams for the simulation parameter subglottal
pressure are depicted for the acoustic pressure evaluated at the microphone
position. A discrimination between the subglottal pressures
$385\text{\,}\mathrm{Pa}$ and
$\\{$775\text{\,}\mathrm{Pa}$,$1500\text{\,}\mathrm{Pa}$\\}$ is possible using
HNR, SPL, Alpha Ratio $\rho_{\alpha}$, and HBI $\eta$, as depicted in Fig. 8
(C), (D), (E), and (F), respectively. Spectral slope allows distinguishing
configurations with $775\text{\,}\mathrm{Pa}$ from
$\\{$385\text{\,}\mathrm{Pa}$,$1500\text{\,}\mathrm{Pa}$\\}$, as it can be
seen from Fig. 8 (A).

Figure 8: Boxplot diagrams of the features spectral slope $a$, CPP
$c_{$2\text{\,}\mathrm{kHz}$}$, HNR $h_{$2\text{\,}\mathrm{kHz}$}$, SPL
$L_{P,$5\text{\,}\mathrm{kHz}$}$, alpha ratio $\rho_{\alpha}$, and HBI $\eta$
as a function of subglottal pressure for acoustic pressure signals at the
microphone position.

####  3.1.1 Intermediate conclusions boxplots

From the boxplot diagrams, statements regarding the significance of the
features’ discrimination ability can be drawn. A good discriminative ability
is provided by using CPP or HNR to discriminate the GC type. The decreasing
trend of CPP for increasing glottal insufficiencies reported in Fig. 6 (B) is
consistent with [57, 46, 73, 74]. In accordance with [24, Fig. 12], where a
distinction of VF motion asymmetries was not evident based on their evaluated
features, the acoustic features of the present work also do not allow for the
discrimination of VF motion asymmetry. Similar to [74, 24], we conclude that
asymmetries in the VF motion are not necessarily related with a reduction in
voice quality. The numerical range of our results is similar to the results of
[24]. Furthermore, the subglottal pressure can be discriminated into two
groups, i.e. $385\text{\,}\mathrm{Pa}$ and
$\\{775,1500\\}$\text{\,}\mathrm{Pa}$$ by using the features SPL, alpha ratio
and HBI, mainly.

###  3.2 LDA/SVM Classification

In the following figures, the background is colored according to the class
boundaries evaluated with SVM. The coloring of the data points indicates if
the SVM is able to classify the point correctly: the inner dot represents the
true class, and the outer ring represents the predicted class by SVM. From
these figures, we see how the simulation configurations cluster in the
dimension-reduced feature space. Furthermore, a graphical indication of class
boundaries is provided. The class boundaries are obtained by classifying
points of a mesh grid in the two-dimensional space with the trained SVM
classifier.

##### LDA & SVM for Glottal Closure Type

In Fig. 9, the result of the dimensionality reduction with LDA is visualized
for all configurations using the pressure signal evaluated at the microphone
point with respect to glottal closure type.

From Fig. 9, we see that a clear separation is possible between GC4 and
$\\{\mathrm{GC1},\mathrm{GC2},\mathrm{GC3}\\}$. The classes GC1, GC2, and GC3
are interwoven in the 2D-space. Thus, the training score222The training score
is obtained as follows. First the SVM is trained with all 24 configurations.
Then the SVM is evaluated at the data points and the percentage of true
classifications in relation to the total data is the training score.  is only
$75\text{\,}\mathrm{\char 37\relax}$. Misclassifications occur primarily for
GC2 and GC3.

Figure 9: (color online) LDA-Reduced dimensionality of feature space and
2D-SVM classification of the glottal closure type for acoustic pressure
signals at the microphone position. The points are colored according to the
true GC type. The colored ring around each point denotes the predicted GC
type. The background colors show the predicted GC types of the dimension-
reduced feature space as labeled, and thus the class boundaries.

##### LDA & SVM for Subglottal Pressure

In Fig. 10, the result of the dimensionality reduction with LDA is visualized
for all configurations using the pressure signal evaluated at the microphone
point with respect to subglottal pressure. In Fig. 10, it is visible that a
good class separability is given for the three subglottal pressures, which is
resembled in the training score of $91.7\text{\,}\mathrm{\char 37\relax}$.
Especially the configurations with $775\text{\,}\mathrm{Pa}$ are well
separated from the others.

Figure 10: LDA-Reduced dimensionality of feature space and 2D-SVM
classification of the subglottal pressure for acoustic pressure signals at the
microphone position. The points are colored according to the true subglottal
pressure. The colored ring around each point denotes the predicted subglottal
pressure. The background colors show the predicted subglottal pressure class
of the dimension-reduced feature space as labeled, and thus the class
boundaries.

##### LDA & SVM for Symmetry

In Fig. 11, the result of the dimensionality reduction with LDA is visualized
for the pressure signal evaluated at the microphone point the microphone
position with respect to the symmetry property. From Fig. 11 we might conclude
that a good class separability is given because the class boundary is clearly
between two clusters. Still, there are many data samples located at the wrong
side of the class boundary which is obvious from the training score of
$70.8\text{\,}\mathrm{\char 37\relax}$.

Figure 11: LDA-Reduced dimensionality of feature space and 2D-SVM
classification of the symmetry for acoustic pressure signals at the microphone
position. The points are colored according to the true symmetry class. The
colored ring around each point denotes the predicted symmetry class. The
background colors show the predicted symmetry class of the dimension-reduced
feature space as labeled, and thus the class boundary.

##  4 Conclusion

Individuals affected by voice disorders have a significant reduction of their
quality of life. Therefore, a simulation model simVoice was developed to
investigate typical characteristics of voice such as sub-glottal pressure and
of functional voice disorders as glottal closure insufficiency and left-right
asymmetry. Based on 24 different voice configurations, we study voice
parameters (HNR, CPP, …) in the prospect of analyzing the flow-induced noise
generation. Although multiple studies are based on experimental data, this is
one of the first studies learning simulated data of healthy and disordered
voice signals. A major limitation of our investigations is that the dataset is
small, but the definition of relevant characteristics are precise and based on
the simVoice model. Furthermore, the small dataset is analyzed by correlation
analysis, and an SVM classifier with RBF kernel. With the use of LDA the
dimensions of the individual representations are visualized in 2D. This allows
to draw first correlations and determine the most important features evaluated
from the acoustic signals in front of the mouth. Based on the definition of
the acoustic features and with focus on the two lower formants, we found that
a low-pass filtering with a cut-off frequency of $2\text{\,}\mathrm{kHz}$
provides a reasonable basis for further analysis. Regarding the results,
conclusions can be drawn as follows. The GC type can be best discriminated
based on CPP. The SPL allows discriminating subglottal pressure into two
groups. The glottal symmetry does not have any effect on the evaluated
features. Using the LDA-dimensionality-reduced feature space, one can best
classify subglottal pressure with 91.7% accuracy. In addition to that the GC
type can be classified with 75% accuracy. Finally, the symmetry property
exhibits interwoven clusters and hence a classification is rather inaccurate.

## Acknowledgment

We acknowledge support from the Austrian Academy of Sciences (ÖAW), research
grant "Understanding voice disorders", received from "Dr. Anton Oelzelt-
Newin’sche Stiftung".

## References

  * [1] I. R. Titze, Principles of Voice Production, 2nd Edition, National Center for Voice and Speech, 2000.
  * [2] S. Schoder, M. Weitz, P. Maurerlehner, A. Hauser, S. Falk, S. Kniesburges, M. Döllinger, M. Kaltenbacher, Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation, J. Acoust. Soc. Am. 147 (2) (2020) 1179–1194.  doi:10.1121/10.0000785.
  * [3] J. B. Park, L. Mongeau, Experimental investigation of the influence of a posterior gap on glottal flow and sound, J. Acoust. Soc. Am. 124 (2) (2008) 1171–1179.  arXiv:https://doi.org/10.1121/1.2945116, doi:10.1121/1.2945116.
  * [4] A. Yamauchi, H. Yokonishi, H. Imagawa, K.-I. Sakakibara, T. Nito, N. Tayama, T. Yamasoba, Age- and gender-related difference of vocal fold vibration and glottal configuration in normal speakers: Analysis with glottal area waveform, J Voice 28 (5) (2014) 525–531.  doi:10.1016/j.jvoice.2014.01.016.
  * [5] M. R. Hoffman, K. Surender, E. E. Devine, J. J. Jiang, Classification of glottic insufficiency and tension asymmetry using a multilayer perceptron, Laryngoscope 122 (12) (2012) 2773–2780.  doi:10.1002/lary.23549.
  * [6] L. A. Rammage, R. C. Peppard, D. M. Bless, Aerodynamic, laryngoscopic, and perceptual-acoustic characteristics in dysphonic females with posterior glottal chinks: A retrospective study, J Voice 6 (1) (1992) 64–78.  doi:10.1016/S0892-1997(05)80010-4.
  * [7] E. C. Inwald, M. Döllinger, M. Schuster, U. Eysholdt, C. Bohr, Multiparametric analysis of vocal fold vibrations in healthy and disordered voices in high-speed imaging, J Voice 25 (5) (2011) 576–590.  doi:10.1016/j.jvoice.2010.04.004.
  * [8] R. R. Patel, A. Dixon, A. Richmond, K. D. Donohue, Pediatric high speed digital imaging of vocal fold vibration: A normative pilot study of glottal closure and phase closure characteristics, Int. J. Pediatr. Otorhinolaryngol. 76 (7) (2012) 954–959.  doi:10.1016/j.ijporl.2012.03.004.
  * [9] B. Schneider, W. Bigenzahn, Influence of glottal closure configuration on vocal efficacy in young normal-speaking women, J Voice 17 (4) (2003) 468–480.  doi:10.1067/S0892-1997(03)00065-1.
  * [10] M. Vaca, I. Cobeta, E. Mora, P. Reyes, Clinical assessment of glottal insufficiency in age-related dysphonia, J Voice 31 (1) (2017) 128.e1–128.e5.  doi:10.1016/j.jvoice.2015.12.010.
  * [11] M. Södersten, S. Hertegård, B. Hammarberg, Glottal closure, transglottal airflow, and voice quality in healthy middle-aged women, J Voice 9 (2) (1995) 182–197.  doi:10.1016/S0892-1997(05)80252-8.
  * [12] S. Moccia, G. O. Vanone, E. De Momi, A. Laborai, L. Guastini, G. Peretti, L. S. Mattos, Learning-based classification of informative laryngoscopic frames, Computer methods and programs in biomedicine 158 (2018) 21–30.
  * [13] H. Cordeiro, J. Fonseca, I. Guimarães, C. Meneses, Hierarchical classification and system combination for automatically identifying physiological and neuromuscular laryngeal pathologies, Journal of voice 31 (3) (2017) 384–e9.
  * [14] M.-H. Laves, J. Bicker, L. A. Kahrs, T. Ortmaier, A dataset of laryngeal endoscopic images with comparative study on convolution neural network-based semantic segmentation, International journal of computer assisted radiology and surgery 14 (3) (2019) 483–492.
  * [15] D. E. Callan, R. D. Kent, N. Roy, S. M. Tasko, Self-organizing map for the classification of normal and disordered female voices, Journal of Speech, Language, and Hearing Research 42 (2) (1999) 355–366.
  * [16] S. N. Awan, N. Roy, Acoustic prediction of voice type in women with functional dysphonia, Journal of voice 19 (2) (2005) 268–282.
  * [17] D. Voigt, M. Döllinger, T. Braunschweig, A. Yang, U. Eysholdt, J. Lohscheller, Classification of functional voice disorders based on phonovibrograms, Artificial Intelligence in Medicine 49 (1) (2010) 51–59.
  * [18] D. Panek, A. Skalski, J. Gajda, R. Tadeusiewicz, Acoustic analysis assessment in speech pathology detection, International Journal of Applied Mathematics and Computer Science 25 (3) (2015) 631–643.
  * [19] S. Umapathy, S. Rachel, R. Thulasi, Automated speech signal analysis based on feature extraction and classification of spasmodic dysphonia: a performance comparison of different classifiers, International Journal of Speech Technology 21 (1) (2018) 9–18.
  * [20] A. Sama, P. Carding, S. Price, P. Kelly, J. Wilson, The clinical features of functional dysphonia, The Laryngoscope 111 (3) (2001) 458–463.
  * [21] P. Schlegel, M. Stingl, M. Kunduk, S. Kniesburges, C. Bohr, M. Döllinger, Dependencies and ill-designed parameters within high-speed videoendoscopy and acoustic signal analysis, Journal of Voice 33 (5) (2019) 811–e1.
  * [22] P. Schlegel, M. Kunduk, M. Stingl, M. Semmler, M. Döllinger, C. Bohr, A. Schützenberger, Influence of spatial camera resolution in high-speed videoendoscopy on laryngeal parameters, PloS one 14 (4) (2019) e0215168.
  * [23] P. Schlegel, S. Kniesburges, S. Dürr, A. Schützenberger, M. Döllinger, Machine learning based identification of relevant parameters for functional voice disorders derived from endoscopic high-speed recordings, scientific reports 10 (1) (2020) 1–14.
  * [24] S. Falk, S. Kniesburges, S. Schoder, B. Jakubaß, P. Maurerlehner, M. Echternach, M. Kaltenbacher, M. Döllinger, 3D-FV-FE aeroacoustic larynx model for investigation of functional based voice disorders, Front. Physiol. 12 (2021) 226.  doi:10.3389/fphys.2021.616985.
  * [25] S. Schoder, P. Maurerlehner, A. Wurzinger, A. Hauser, S. Falk, S. Kniesburges, M. Döllinger, M. Kaltenbacher, Aeroacoustic sound source characterization of the human voice production-perturbed convective wave equation, Appl. Sci. 11 (6).  doi:10.3390/app11062614.
  * [26] H. Sadeghi, M. Döllinger, M. Kaltenbacher, S. Kniesburges, Aerodynamic impact of the ventricular folds in computational larynx models, J. Acoust. Soc. Am. 145 (4) (2019) 2376–2387.  doi:10.1121/1.5098775.
  * [27] H. Sadeghi, S. Kniesburges, M. Kaltenbacher, A. Schützenberger, M. Döllinger, Computational models of laryngeal aerodynamics: Potentials and numerical costs, J Voice 33 (4) (2019) 385–400.  doi:10.1016/j.jvoice.2018.01.001.
  * [28] H. Sadeghi, S. Kniesburges, S. Falk, M. Kaltenbacher, A. Schützenberger, M. Döllinger, Towards a clinically applicable computational larynx model, Appl. Sci. 9 (11).  doi:10.3390/app9112288.
  * [29] P. Maurerlehner, S. Schoder, C. Freidhager, A. Wurzinger, A. Hauser, F. Kraxberger, S. Falk, S. Kniesburges, M. Echternach, M. Döllinger, M. Kaltenbacher, Efficient numerical simulation of the human voice, e & i Elektrotechnik und Informationstechnik 138 (3) (2021) 219–228.  doi:10.1007/s00502-021-00886-1.
  * [30] F. Alias, J. C. Socoro, X. Sevillano, A review of physical and perceptual feature extraction techniques for speech, music and environmental sounds, Appl. Sci. 6 (5).  doi:10.3390/app6050143.
  * [31] F. Eyben, Real-time Speech and Music Classification by Large Audio Feature Space Extraction, Springer Theses, Springer International, 2016.  doi:10.1007/978-3-319-27299-3.
  * [32] J. T. Geiger, B. Schuller, G. Rigoll, Large-scale audio feature extraction and svm for acoustic scene classification, in: 2013 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 2013, pp. 1–4.  doi:10.1109/WASPAA.2013.6701857.
  * [33] I. Mierswa, K. Morik, Automatic feature extraction for classifying audio data, Mach. Learn. 58 (2) (2005) 127–149.  doi:10.1007/s10994-005-5824-7.
  * [34] S. Schoder, F. Kraxberger, S. Falk, A. Wurzinger, K. Roppert, S. Kniesburges, M. Döllinger, M. Kaltenbacher, Error detection and filtering of incompressible flow simulations for aeroacoustic predictions of human voice, J. Acoust. Soc. Am. (under review).
  * [35] M. Kaltenbacher, Numerical Simulation of Mechatronic Sensors and Actuators: Finite Elements for Computational Multiphysics, 3rd Edition, Springer, Berlin–Heidelberg, 2015.  doi:10.1007/978-3-642-40170-1.
  * [36] S. Schoder, K. Roppert, opencfs: Open source finite element software for coupled field simulation–part acoustics, arXiv preprint arXiv:2207.04443.
  * [37] B. Kaltenbacher, M. Kaltenbacher, I. Sim, A modified and stable version of a perfectly matched layer technique for the 3-d second order wave equation in time domain with an application to aeroacoustics, J. Comput. Phys. 235 (2013) 407–422.  doi:10.1016/j.jcp.2012.10.016.
  * [38] J. Probst, A. Lodermeyer, S. Fattoum, S. Becker, M. Echternach, B. Richter, M. Döllinger, S. Kniesburges, Acoustic and aerodynamic coupling during phonation in mri-based vocal tract replicas, Appl. Sci. 9 (17).  doi:10.3390/app9173562.
  * [39] S. Falk, Numerical computation of functional voice disorders, Ph.D. thesis, FAU Erlangen-Nürnberg (2021).
  * [40] M. A. Belsky, S. D. Rothenberger, A. I. Gillespie, J. L. Gartner-Schmidt, Do phonatory aerodynamic and acoustic measures in connected speech differ between vocally healthy adults and patients diagnosed with muscle tension dysphonia?, J Voice 35 (4) (2021) 663.e1–663.e7.  doi:10.1016/j.jvoice.2019.12.019.
  * [41] O. Garaycochea, J. M. A. Navarrete, B. del Río, S. Fernández, Muscle tension dysphonia: Which laryngoscopic features can we rely on for diagnosis?, J Voice 33 (5) (2019) 812.e15–812.e18.  doi:10.1016/j.jvoice.2018.04.015.
  * [42] S. Kniesburges, A. Lodermeyer, M. Semmler, Y. K. Schulz, A. Schützenberger, S. Becker, Analysis of the tonal sound generation during phonation with and without glottis closure, J. Acoust. Soc. Am. 147 (5) (2020) 3285–3293.  doi:10.1121/10.0001184.
  * [43] R. Patel, D. Dubrovskiy, M. Döllinger, Characterizing vibratory kinematics in children and adults with high-speed digital imaging, J Speech Lang Hear Res 57 (2) (2014) S674–S686.  doi:10.1044/2014_JSLHR-S-12-0278.
  * [44] Z. Zhang, Compensation strategies in voice production with glottal insufficiency, J Voice 33 (1) (2019) 96–102.  doi:10.1016/j.jvoice.2017.10.002.
  * [45] A. Lodermeyer, S. Becker, M. Döllinger, S. Kniesburges, Phase-locked flow field analysis in a synthetic human larynx model, Exp. Fluids 56 (4) (2015) 77\.  doi:10.1007/s00348-015-1942-6.
  * [46] V. Birk, S. Kniesburges, M. Semmler, D. A. Berry, C. Bohr, M. Döllinger, A. Schützenberger, Influence of glottal closure on the phonatory process in ex vivo porcine larynges, J. Acoust. Soc. Am. 142 (4) (2017) 2197–2207.  doi:10.1121/1.5007952.
  * [47] P. Maurerlehner, Efficient FEM model of human phonation, Master’s thesis, TU Wien, Vienna (Mar. 2020).
  * [48] P. Woo, A. K. Parasher, T. Isseroff, A. Richards, M. Sivak, Analysis of laryngoscopic features in patients with unilateral vocal fold paresis, Laryngoscope 126 (8) (2016) 1831–1836.  doi:10.1002/lary.25790.
  * [49] S. S. Azar, P. Pillutla, L. K. Evans, Z. Zhang, J. Kreiman, D. K. Chhetri, Perceptual evaluation of vocal fold vibratory asymmetry, Laryngoscope 131 (12) (2021) 2740–2746.  doi:10.1002/lary.29679.
  * [50] A. V. Oppenheim, J. R. Buck, R. W. Schafer, Discrete-Time Signal Processing, 3rd Edition, Prentice-Hall signal processing series, Pearson Education, 2014.  doi:10.5555/1795494.
  * [51] E. Yumoto, W. J. Gould, T. Baer, Harmonics-to-noise ratio as an index of the degree of hoarseness, J. Acoust. Soc. Am. 71 (6) (1982) 1544–1550.  doi:10.1121/1.387808.
  * [52] K. Lee, C. Moon, Y. Nam, Diagnosing vocal disorders using cobweb clustering of the jitter, shimmer, and harmonics-to-noise ratio, KSII T. Internet Info. 12 (11) (2018) 5541–5554.  doi:10.3837/tiis.2018.11.020.
  * [53] J. Kreiman, B. R. Gerratt, M. Garellek, R. Samlan, Z. Zhang, Toward a unified theory of voice production and perception, Loquens 1 (27135054) (2014) e009.  doi:10.3989/loquens.2014.009.
  * [54] C. T. Ferrand, Harmonics-to-noise ratio: an index of vocal aging, J Voice 16 (2002) 480–7.  doi:10.1016/s0892-1997(02)00123-6.
  * [55] B. P. Bogert, M. J. R. Healy, J. W. Tukey, The quefrency alanysis of time series for echoes: Cepstrum, pseudo-autocovariance, cross-cepstrum and saphe cracking, in: M. Rosenblatt (Ed.), Proc. of the Symposium on Time Series Analysis, John Wiley & Sons, Inc., New York, 1963, pp. 209–243.
  * [56] J. Hillenbrand, R. A. Cleveland, R. L. Erickson, Acoustic correlates of breathy vocal quality, J. Speech Hear. 37 (1994) 769–778.  doi:10.1044/jshr.3704.769.
  * [57] J. Hillenbrand, R. A. Houde, Acoustic correlates of breathy vocal quality: dysphonic voices and continuous speech., J. Speech Hear. 39 (1996) 311–321.  doi:10.1044/jshr.3902.311.
  * [58] P. Welch, The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, IEEE Trans. Audio Electroacoust. 15 (2) (1967) 70–73.  doi:10.1109/TAU.1967.1161901.
  * [59] A. L. Rosenthal, S. Y. Lowell, R. H. Colton, Aerodynamic and acoustic features of vocal effort, J Voice 28 (2) (2014) 144–153.  doi:10.1016/j.jvoice.2013.09.007.
  * [60] Y. Maryn, P. Corthals, P. Van Cauwenberge, N. Roy, M. De Bodt, Toward improved ecological validity in the acoustic measurement of overall voice quality: Combining continuous speech and sustained vowels, J Voice 24 (5) (2010) 540–555.  doi:10.1016/j.jvoice.2008.12.014.
  * [61] Y. Maryn, M. De Bodt, N. Roy, The acoustic voice quality index: Toward improved treatment outcomes assessment in voice disorders, J. Commun. Disord. 43 (3) (2010) 161–174.  doi:10.1016/j.jcomdis.2009.12.004.
  * [62] Y. Maryn, N. Roy, M. De Bodt, P. Van Cauwenberge, P. Corthals, Acoustic measurement of overall voice quality: a meta-analysis., J. Acoust. Soc. Am. 126 (2009) 2619–34.  doi:10.1121/1.3224706.
  * [63] A. Behrman, D. Finan, Speech and voice science, 4th Edition, Plural Publishing, 2021\.
  * [64] D. P. W. Ellis, An Introduction to Signal Processing for Speech, 2nd Edition, Blackwell Handbooks in Linguistics, Wiley-Blackwell, 2010, pp. 757–780.  doi:10.1002/9781444317251.ch20.
  * [65] B. Gold, N. Morgan, D. Ellis, Speech and Audio Signal Processing, 2nd Edition, John Wiley & Sons, 2011.  doi:10.1002/9781118142882.
  * [66] J. Kreiman, Y.-L. Shue, G. Chen, M. Iseli, B. R. Gerratt, J. Neubauer, A. Alwan, Variability in the relationships among voice quality, harmonic amplitudes, open quotient, and glottal area waveform shape in sustained phonation, J. Acoust. Soc. Am. 132 (23039455) (2012) 2625–2632.  doi:10.1121/1.4747007.
  * [67] G. Chen, J. Kreiman, Y. Shue, A. Alwan, Acoustic correlates of glottal gaps, in: Proc. INTERSPEECH 2011, ISCA, Florence, Italy, 2011, pp. 2673–2676.
  * [68] S. Kakouros, O. Räsänen, P. Alku, Comparison of spectral tilt measures for sentence prominence in speech—effects of dimensionality and adverse noise conditions, Speech Commun. 103 (2018) 11–26.  doi:10.1016/j.specom.2018.08.002.
  * [69] B. Hammarberg, B. Fritzell, J. Gaufin, J. Sundberg, L. Wedin, Perceptual and acoustic correlates of abnormal voice qualities, Acta Otolaryngol. 90 (1-6) (1980) 441–451.  doi:10.3109/00016488009131746.
  * [70] S. Patel, K. Scherer, J. Sundberg, E. Björkner, Acoustic markers of emotions based on voice physiology, in: Proc. Speech Prosody, Vol. 100865, ISCA, 2010.
  * [71] C. M. Bishop, Pattern Recognition and Machine Learning, Information Science and Statistics, Springer, 2006.  doi:10.5555/1162264.
  * [72] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, J Mach Learn Res 12 (2011) 2825–2830.  doi:10.5555/1953048.2078195.
  * [73] F. Thornton, M. Döllinger, S. Kniesburges, D. Berry, C. Alexiou, A. Schützenberger, Impact of subharmonic and aperiodic laryngeal dynamics on the phonatory process analyzed in ex vivo rabbit models, Appl. Sci. 9 (9).  doi:10.3390/app9091963.
  * [74] Z. Zhang, J. Kreiman, B. R. Gerratt, M. Garellek, Acoustic and perceptual effects of changes in body layer stiffness in symmetric and asymmetric vocal fold models, J. Acoust. Soc. Am. 133 (1) (2013) 453–462.  doi:10.1121/1.4770235.