Variability and Central Tendencies in Speech Production

Speech is notoriously variable, but our understanding of this variability continues to evolve. Variability has typically been taken as an indication of failure to reach a desired target due to physical or neurological limits. However, it is likely that some variability is beneficial, an effect that has been found in other domains. Part of the effort to separate beneficial from destructive variability must be to understand the distribution of values around a speech target. One aspect that is commonly measured is the standard deviation of some objective aspect of speech. The standard deviation is most meaningful for normal distributions, and the assumption in speech research has been that values are indeed normally distributed. This has not been rigorously tested, however, as the test of normality requires a large number of samples (some studies suggest a minimum of 200) to determine whether the data is normally distributed or not. Speech research (and, indeed, most research with humans) seldom reaches such numbers for a consistent environment. Here, an initial estimate for 300 repetitions of English words by a single speaker are presented. The words were pseudo-randomized with an equal number of filler items, so that immediate repetitions (and the neural and physical fatigue repetition can cause) were avoided. 100 trials were collected on each of three days. Words were chosen to have very little coarticulatory influence (“heed,” “ode”/“owed”) or sizable coarticulatory influence (“geek,” “dote”). Measurements of vowel formants at acoustic midpoints indicated that the distributions were indeed normal. This was true even of the high coarticulatory environment, which some theories would predict would be skewed by the vowel’s reaching the edge of an acceptable region. The current results indicate that vowel targets are consistent for different environments. Further, the range of the distributions was quite similar across the two types of environment, being, for example, about 100 Hz for F1. The amount of variability is fairly substantial but can be presumed to be beneficial, as all items were heard correctly. The normality of the distribution nonetheless indicates a control structure that accommodates the coarticulatory environment at the level of planning.

Speech is notoriously variable, but our understanding of this variability continues to evolve. Variability has typically been taken as an indication of failure to reach a desired target due to physical or neurological limits. However, it is likely that some variability is beneficial, an effect that has been found in other domains. Part of the effort to separate beneficial from destructive variability must be to understand the distribution of values around a speech target. One aspect that is commonly measured is the standard deviation of some objective aspect of speech. The standard deviation is most meaningful for normal distributions, and the assumption in speech research has been that values are indeed normally distributed. This has not been rigorously tested, however, as the test of normality requires a large number of samples (some studies suggest a minimum of 200) to determine whether the data is normally distributed or not. Speech research (and, indeed, most research with humans) seldom reaches such numbers for a consistent environment. Here, an initial estimate for 300 repetitions of English words by a single speaker are presented. The words were pseudo-randomized with an equal number of filler items, so that immediate repetitions (and the neural and physical fatigue repetition can cause) were avoided. One hundred trials were collected on each of 3 days. Words were chosen to have very little coarticulatory influence ("heed," "ode"/"owed") or sizable coarticulatory influence ("geek," "dote"). Measurements of vowel formants at acoustic midpoints indicated that the distributions were indeed normal. This was true even of the high coarticulatory environment, which some theories would predict would be skewed by the vowel's reaching the edge of an acceptable region. The current results indicate that vowel targets are consistent for different environments. Further, the range of the distributions was quite similar across the two types of environment, being, for example, about 100 Hz for F1. The amount of variability is fairly substantial but can be presumed to be beneficial, as all items were heard correctly. The normality of the distribution nonetheless indicates a control structure that accommodates the coarticulatory environment at the level of planning.

INTRODUCTION
Variability is a well-known feature of speech, as it is with other biological systems. Although excessive variability can signify lack of motor control, lack of variability can itself be pathological (e.g., Dinstein et al., 2015). Variability in input has been shown to be helpful in learning (e.g., Bradlow et al., 1997;Preston et al., 2018), and variability in production can give a range of options for adapting to novel situations (e.g., Ossmy et al., 2018).
A typical assumption is that variability is normally distributed; indeed, it is typical to the point that most studies do not explicitly state that assumption. The successful analysis of results in those studies suggests that the assumption is justified to a great extent. Although many statistical tests provide replicable results even when their requirements are violated (Lix et al., 1996), there are indications that other results can be greatly affected (Cain et al., 2017). Concerns over the effects of skew and kurtosis have motivated the move to linear mixed effects models (Baayen et al., 2008;Pouplier et al., 2017), but the meaning of the distributions themselves is not addressed by such analyses.
Non-normal distributions of data can indicate that more than one process is affecting the distribution [A normal distribution can arise from multiple sources if the samples are independent and identically distributed (central limit theorem)]. If the distribution is bimodal, then the single measurement under consideration may be treating two effects as if they were one. If there is skew in the distribution, there may be influences at work that need to be addressed. For statistical purposes, skew may invalidate some tests, such as ANOVA (e.g., Harwell et al., 1992). The more interesting effect is that it may indicate an influence on the behavior under study. The most common of these, of course, is a boundary effect, when the mean of a distribution is close to a physical limit, that is, when the standard deviation simply cannot extend as far as it would without the constraint. Both of these effects can be informative rather than a hindrance to analysis if they are examined on their own. That is one purpose of the present experiment.
In this study, we examined the pattern of articulatory variability in vowel targets for English. Although direct articulatory measurements are more readily obtained now than in previous years, they are still demanding in data collection and analysis, making them challenging for large-scale studies [though see discussion below of a physiological study in Tilsen (2017)]. Here, we needed many repetitions in order to examine the distributional characteristics of the productions, and so we, like many others, relied on the acoustic output to index the articulatory activity. Not only is the acoustic output reliably shaped by the articulation (Fant, 1960;Iskarous, 2010), there is also evidence that variability in the acoustic domain is highly related to the variability in the articulatory domain (Whalen et al., 2018). The use of acoustics therefore is a reasonable first step in analyzing production variability.
The focus here is on random variability, not structured variability, so we needed to focus on single targets. There is a great deal of structured variability due to vocal tract length, coarticulation, emotion, etc. (e.g., Best, 2015), and such variation is of great importance for understanding the entirety of the phonetic system. If it were possible to code all of those structured effects on formant values, we might be able to assess distributions from large speech corpora; the residual after removing the structured effects would be the unstructured variability. However, no corpora are annotated to that extent, and it may be that none ever will be. The number of systematic sources of variability is sizable and generally expanding as more studies are completed. Relying on our accurate account of those factors is not possible at present, given, for example, the relatively inaccurate methods for vocal tract normalization (e.g., Flynn, 2011). Thus, we relied on multiple repetitions of non-contextualized words by a single speaker.
The level of variability due to motor noise and other intrinsic factors must be examined with productions that lack, to the extent possible, structured variability. "Intrinsic" factors are here conceptualized as distinct from structured ones, and they would include such variables as arousal state, location in the breath cycle, and changes in the motor program (either "intentional" or not). The boundary between intrinsic and systematic is not firm, however, and they may not really affect variability differently. We nonetheless wanted to avoid as many factors unrelated to the motor program as we could. To that end, we elicited multiple repetitions of target words so that, ideally, only variability in motor planning and execution remained. Fatigue of motor systems in sustained repetition is well-attested even if the underlying cause (central nervous system, the neuromuscular junction, or metabolic changes in the muscle fiber) is difficult to ascertain (e.g., Bigland-Ritchie, 1981). Thus, the paradigm of having a speaker produce many repetitions of a word [such as the 1,000 sequential repetitions of the word "bucket" in Kello et al. (2008)] can be expected to induce variability based on sheer physical and neural fatigue that are not relevant to understanding what speakers do when they are producing their ideal version of a word. We therefore collected our target words in lists which contained an equal number of filler items, allowing the neurons and muscles to reset between productions.
Direct instructions to eliminate variability do not appear to be successful and may even be counterproductive. In a study of multiple repetitions of target items, Tilsen (2017) provided feedback about consistency in an attempt to eliminate variability. It failed: Speakers continued to have variability, and the variability was structured across independent motor systems. For our purposes, the results indicated that providing feedback about individual productions was not effective in eliminating variability and therefore increased the cognitive load on the speaker without necessarily modifying the speaker's behavior. We therefore strove for consistency simply by asking the speaker to be consistent.
Formant measurements are known to be influenced by fundamental frequency (F0), but large datasets require automatic measurements that currently include such influences. Vowels are well-described by the formants (Fant, 1960), but it is really the resonances that are the true object of interest (Titze et al., 2015). Acoustic formant analysis tends to follow the most intense harmonic near a resonance (F0-effect) (Klatt, 1986;Shadle et al., 2016), but listeners respond to the true resonances, not the measured formants (Klatt, 1986). In the present study, we found that F0 effects were minimal due to the great consistency of F0 by our speaker, so that automatic measurements of formants were usable.
Many tokens are required to analyze the distribution of variability, but studies of speech seldom obtain the required amount. If 20 tokens are collected, we can obtain a fairly defensible estimate of the central tendency (mean) of the distribution, but a sample of only 20 tokens will almost always appear to be normally distributed, even if the true distribution is not normal. Mardia (1970) found that there was more than twice as much evidence for either atypical skewness or kurtosis when the sample size exceeded 106 (46 vs. 94%), indicating that large samples are needed for these measures. In a simulation of models with many parameters, Lerche et al. (2017) found that 200 trials provided good estimations for three-and four-parameter models (p. 522). These are not exact matches to the current experiment, but they give an indication of how many trials can be expected to give solid results. Thus, a sample size of 200 should provide good evidence of distributional properties; we oversampled by obtaining 300 repetitions.
Two environments were studied, allowing us to study intrinsic variability in two extrinsic changes in coarticulation. The first was an /(h)Vd/ environment, which has been shown to have small if any effects on vowel midpoint formant measures in comparison to isolated vowels (Stevens and House, 1963;Ohde and Sharf, 1975). The second was an environment of consonants that differed maximally from the vowel's position, that is Our first analysis contrasts two hypotheses about the effect of coarticulation on the distribution of formant values. The first hypothesis, based on the "window" model of coarticulation (Keating, 1990) is that a neutral environment would have small skewness values while a coarticulated one would have larger skewness. The alternative model, labeled more generically as the non-window model, predicts non-skewed distributions for both environments. The rationale for this can be seen in Figure 1. The window model hypothesizes that the planning stage contains no central target for a segment, only a range (the window) of variability. The implementation is then the result of an interpolation process that finds an optimal path through connected windows with minimum articulatory effort. A window is defined as a pair of minimum and maximum values in a physical dimension that the observed productions are bounded by Keating (1990, pp. 455-456). Thus, a boundary effect on the skewness of distribution should occur if the path from one window to another is most easily accomplished by moving close to an edge. The predictions of Guenther's (1995) "convex region theory" would seem to be the same as the window model's, because the region is meant to be sufficient for the production of a target. His regions are multidimensional and include somatosensory space, so the acoustic predictions are not straightforward. Nonetheless, because the theory is meant to account for such features as vowel reduction (undershoot) (Lindblom, 1963(Lindblom, , 1983, it would seem that it would make the same prediction as the window model in this case: Vowels must enter the convex region to be successful, so there should be few productions outside the convex region. Productions that enter the region more deeply will also be successful, but less common. Formant values would therefore be expected to show a skewed distribution. A further complication is that segments have somatosensory targets as well as acoustic ones, resulting in separate error calculations for each (e.g., Terband et al., 2009). Whether this later interaction would affect the distribution has not been tested. Figure 1 shows, schematically for F2 alone, the executed path of F2 (red solid lines) for /o w d/ and /do w t/, necessarily the same in both the window model and the opposing non-window model. Hypothetical resultant F2 trajectories are shown from the onset of the syllable to the end of the first component of the vowel (omitting the offglide /w/ and the stop coda). The difference in the models is the control parameters, shown by the black dotted lines (the range of target). For the window model, this range defines planning parameters that are the same regardless of context. For the non-window model, this range represents a confidence interval of a normal distribution generated by, for example, a non-linear dynamical system (Saltzman and Munhall, 1989). For /o w d/ (Figure 1A), both the window and non-window models predict that the F2 trajectory will concentrate in the middle of the target range without skewness. For /do w t/, the window model predicts that the F2 trajectory in the onset will accommodate the desired minimum effort by being toward the lower boundary of the target range, resulting in positive skew. In the middle of /o/, the skewness will be negative as the path enters the upper part of the range (Figure 1B). The non-window model predicts that the F2 distributions for /do w t/ should be normal throughout the whole trajectory ( Figure 1C). Note that the two predictions have the same central tendency of the output trajectory but different predicted patterns of skewness. This is because the target range for a segment in the window model is always the same in all contexts, while in the non-window model the target range can be variable in different contexts as the result of gestural interactions. We chose /o/, despite its known diphthongal offglide (Pike, 1947), because the mid vowels have less chance of abutting a physiological limit, as we expect for /i/.
The second hypothesis is that /i/ should exhibit formant distributions that are somewhat skewed (i.e., positive skew in F1), given that the constriction for /i/ is limited as it approaches the hard palate.
Formants for vowels are rarely stable throughout the vocalic segment, whether the vowel is perceived as diphthongal or not (Hillenbrand and Nearey, 1999). Our analysis examines both the midpoint of the vocalic segment, often seen at the target of the vowel, and the trajectories as well.

EXPERIMENT
Many repetitions of linguistic utterances are needed to address the issue of the normality of the distributions of vowel formants. This need dictated that the target words be produced in isolation so that the recording sessions would be short enough to be tolerable by the speaker. Filler items were needed to avoid excessive repetition and its concomitant shift in neural and muscular response.

Speaker
The speaker was a native speaker of American English. He is a trained phonetician as well as an instructor for the singing voice. He provided written informed consent as approved by the CUNY University Integrated IRB (City University of New York).

Materials
The target words were "heed, " "geek, " "owed"/"ode, " and "dote." The homophones "owed" and "ode" were used as a condition for the experiment that was addressed by the filler items (not discussed here). Results for those items will be presented both separately and combined. Filler words were 25 homophones such as "air"/"ere" and "plain"/"plane."

Procedure
Recordings were made in a sound-attenuated booth at the Graduate Center of the City University of New York (CUNY). A free-field microphone (PCB Piezotronics 482C16) with builtin pre-amp (PCB Piezotronics 378B02) was used. An AD Instruments Power Lab (8/35-1008) data acquisition device with a Dell Optiplex 9010 computer processed signals, which were sampled at a rate of 44.1 kHz.
Recordings were made on 3 separate days, separated by 4 months in the first case and 18 months in the second. Each target word (or two words, in the case of "ode"/"owed") occurred 100 times in the randomized list for each day. Each group of 8 items contained one example of each target word (with "ode" and "owed" randomly assigned) along with 4 filler items. The 50 filler items were randomized twice, once for the first half of the session and another time for the second half.
Words were presented in standard orthography, one at a time, on a computer screen controlled by the Presentation program (https://www.neurobs.com/).

Measurements
The recorded audio files were downsampled to 16 kHz and forced-aligned via FAVE-align tool (Rosenfelder et al., 2014), then manually corrected when necessary. Formant frequencies were measured by the Burg method of linear predictive coding (LPC) (window size = 45 ms; step size = 2 ms, number of poles = 14, pre-emphasis from 50 Hz; Nyquist frequency = 5,000 Hz) with Viterbi tracking using Praat (Boersma and Weenink, 2019). The tracked formant frequencies were time-normalized into 11 points, representing measurements from 0 to 100% in steps of 10%.

Statistics
We carried out univariate normality tests on static formant values, separately for each word produced in each day. We ran the Shapiro et al. (1968) tests of skewness and kurtosis using the "normtest" package in R, taking F1 and F2 as separate dependent variables. To control for the inflation of Type I error due to multiple hypothesis testing, p-values of normality tests were adjusted by Benjamini and Hochberg's (1995) approach of "False Discovery Rate" (FDR). To predict the dynamic formant patterns from the data, we fit a Smoothing Spline ANOVA (SSANOVA) (Gu, 2002) model by adding "Context" and "Time" factors and their interaction, with a random effect of "Day, " separately for each vowel and each formant.

Vowel Midpoint
In order to make the initial analysis tractable, a single time point was used: 50% of the duration of the vocalic segment. Table 1 summarizes the formant values along with their standard deviations (SDs). Results are shown for the 3 recording days separately as well as for the four forms (/hid/, /gik/, /o w d/, and /do w t/) across days. Values for "ode" and "owed" are combined in the form /o w d/.

Distributions and Normality Tests
Normality of the formant distributions was tested statistically, but it can also be visually represented by the kernel density estimation (KDE). Figure 2 presents the distributions of F1 (left column) and F2 (right column) for the four forms (in each row) separately for each day (Day 1: blue solid lines; Day 2: red dotted lines; Day 3: green dashed lines). As can be seen in Figure 2, the distributions were quite regular for each day (i.e., 100 repetitions of each target form), but each day was somewhat different. From Figure 2, we can observe skewness on F2 distribution for /gik/ and /o w d/ produced in Day 1 and for /hid/ in Day 2, as well as on F1 distribution for /gik/ in Day 3. Table 2 summarizes the statistics of the moment coefficient of skewness and the excess kurtosis based on the distributions of F1 and F2 values measured at the vowel midpoint (Figure 2). Excess kurtosis is calculated as kurtosis (the fourth moment) minus three. The expected values for both skewness and excess kurtosis are zero for a normal distribution. Positive values of skewness indicate that the distribution was higher than the mean more often than expected (longer tail in higher frequency). An absolute value of skewness >1 is considered as highly skewed, and an absolute value in between 0.5 and 1 indicates moderately skewed. Positive excess kurtosis indicates the distribution is "skinnier" than a normal distribution with "fatter" tail presumably due to outliers, while negative excess kurtosis indicates the opposite. The indications of significance symbols were based on the p-values adjusted by Benjamini and Hochberg's (1995) FDR method for each block. For example, in the top-left block (F1 for /hid/) of Table 1, the six p-values (not shown) for the tests of both skewness and kurtosis in the 3 days were entered into FDR-adjustment; the (family-wise) null hypothesis is that none of the six statistics came from a normal distribution; any one FDR-adjusted p-value in a block that meets the significance level suggests rejection of such null hypothesis. The statistics in Table 2 showed that the distribution of F1 for /gik/ produced in Day 1 and those of F2 for /gik/ and /o w d/ produced in Day 1 are significantly skewed, which conformed to the shapes of distributions observed in Figure 2.
The dynamic pattern of skewness makes the evidence for an effect on the distributions even less likely. In Figure 3, skewness is calculated for each of the 11 time points of the time-normalized data. Solid circles indicate significant skewness while empty circles non-significant. Significance was based on FDR-adjusted p-values across 11 points of skewness separately for each day and for each formant, with a family-wise null hypothesis as none of the measured values of skewness in the 11 points conforms to normal distribution. Because the consonant(s) at the syllable boundary should have windows of their own, the skew could be expected to change over the course of the syllable, perhaps with a midpoint differing from both ends. Such a pattern is seen for F2 of /o w d/ on day 1. However, two aspects of that pattern are inconsistent with our predictions: The onset of /o w d/ should not be skewed, given that the target can be achieved from the beginning of the utterance. Even if there were an explanation for the presence of the skew, there is no obvious reason that the skew would not be present throughout the vocalic segment (up until the transitions for the final stop). Days 2 and 3, as can be seen, had radically different patterns; there is no clear interpretation for the differences. In short, whatever was skewing some of the formant distributions on some days was not systematic enough to be explained by either the window model or by the non-window model (see Figure 1). Figure 4 further visualizes the distributions of formant frequencies for all time points. Each gray-scaled contour represents the KDE-estimated probability density function (as those distributions displayed in Figure 2) at each time point; darker color indicates higher probability. Red crosses track the means of the distributions along the time course, and blue circles the mode (estimated by measuring the peak of probability density function) of distributions. The difference between mean and mode is known as the nominator of Pearson's mode skewness [(mean-mode)/SD]: If the mean is higher than the mode, it indicates positive skewness, which is a conservative visualization of the direction of skewness. Note that mode skewness may not be perfectly consistent with moment coefficient of skewness (as in Table 2). Figure 4 is largely consistent with Figure 3 and provides more information of probability distributions of formant values at each time point.

Dynamic Formant Patterns
The changes in formant location for the words across all 3 days were examined. The time-normalized values were used. A smoothing spline ANOVA (SSANOVA) was computed separately for F1 and F2 for /hid/ vs. /gik/ (Figure 5) and for /o w d/ vs. /do w t/ (Figure 6). In such displays, the 95% Bayesian confidence intervals (shown in color around the mean formant values) are assumed to be statistically different when they do not overlap. The amount of divergence that is needed before the result is "significant" is debatable, but the existence of a visually distinct region suggests that the trajectories do differ in some ways. As can be seen in Figures 5, 6, the first two formants were constantly changing, leaving no portion that was truly "steadystate." Indeed, inclusion of such minor variability has been shown to improve identification of synthetic versions of the target vowels (Hillenbrand and Nearey, 1999). Other predictable FIGURE 2 | F1 and F2 distributions for the four words (measured at the vowel midpoint) separately for each day. aspects appeared. A separate SSANOVA (not presented here) comparing the homonyms "ode" and "owed" showed that they were, indeed, virtually identical. The formants for the shared alveolar stop at the end of the /o w / words converged (Figure 6).
The formants for the distinct places of articulation of the final stops for the /i/ words diverged (Figure 5). F2 was distinguished at the final portion of the trajectory in Figure 5 and the first half of the trajectory in Figure 6. What was perhaps somewhat surprising was the overall dissimilarity of F1 for the two contexts for the /o w / words but not for the /i/ words. Still, the differences were small (45 Hz for /o w / words, and 2 Hz for /i/ words at the midpoint). Although "geek" was intended to have velar productions on either side of the vowel, the low F2 values at onset indicate that this speaker used a very fronted place of articulation for the initial stop. Thus, the F2 pattern was quite linear, while the F2 of "dote" (Figure 6) behaved as intended. The vowel of "ode"/"owed" was, as expected, rather diphthongal, with F1 changing by about 65 Hz from time points 4 to 8 (the likely limits of coarticulatory effects of the stop). The vowel of "heed, " by contrast, changed by about 10 Hz over those same time points.

DISCUSSION
Multiple repetitions of English words in a fairly isolated state were found to have formants that were only slightly different from normality. Having a sizable number of tokens is necessary for such an analysis, but the biological constraints on speakers make collection challenging. Here, we reduced the constraints as much as possible by interleaving the tokens with filler items, but that limited us to collecting 100 repetitions in any one session.  As can be seen in Figures 2, 3, the formants obtained were quite consistent within those sessions; the small differences across sessions were smaller than the likely measurement error of the LPC analysis. Although changes in articulation across different days or even time of day (Heald and Nusbaum, 2015) have been reported before, the differences here are negligible.
The window model hypothesis that coarticulation would skew the distributions was not supported, while the non-window model was consistent with the lack of skewness. The trajectories were normally distributed not only near the midpoint of the vocalic segment, but throughout the production (Figures 3, 4). Such a result is inconsistent with the "window" model in which the motor plan contains only target regions and not trajectories; in execution, segments reach the edge of one target region before moving on to the next (Keating, 1990). It is consistent with nonwindow model in which a motor plan takes the entire context into account from the beginning; overlapping activations for the gestures or segments then unfold in execution in such a way that  variability is structured by the interactions of the overlapping control parameters of gestures or segments.
The second hypothesis, that the /i/ formants would have skewed distributions because of the boundary effect of the hard palate, was not supported. Not only were there very few individual time points with significant skew, there was no discernable pattern to the skew either. For this speaker, at least, the constraints on articulation of the high front vowel were well-accommodated, so that the distributions of formants were unaffected by the physiological limits. Standard deviations were small but non-negligible (at midpoint, for /i/, 5.8% for F1, 2.3% for F2; for / o w /, 4.6% for F1, 3.8% for F2). It would seem that there is enough variability for a skewed distribution to be evident, if it were present. Instead, the formant distributions appear to be normal through the duration of the syllable.
Future studies are desirable to explore these issues further. Only one speaker was analyzed here, and he was chosen in part for his many years of practice and instruction in broadcast speaking. The resulting consistency was useful for having manageable amounts of variability, but less skilled speakers may show different patterns. Indeed, the kinds of variability that result may differ by such factors as speech sound disorder or speaking in a second language. Other acoustic or articulatory measures could be made, although the strongest predictions in the field have been about formant values. Measuring variability across the vowel system rather than for just two vowels would be useful (Whalen et al., 2018), although the number of tokens required becomes rather large. Finding word tokens that maintain the voicing of the final consonant would also be desirable. Other statistical approaches, such as Generalized Additive Mixed Models, may provide further insight.
Overall, the results for this speaker support the use of statistics that rely on normal distributions for analyzing formant values. As such, the results also support the use of Gaussian priors in Bayesian linear mixed models (Vasishth et al., 2018). Using the results of a single speaker has intrinsic drawbacks, so the current results can only be preliminary. Further, the formant values themselves are subject to many measurement errors (Klatt, 1986;Shadle et al., 2016), but, within those limits, estimation of the central tendencies for formants are relatively good, at least for F0s <200 Hz (Chen et al., 2019). The present data did not support models that assume target regions; instead, entire trajectories were normally distributed throughout the vocalic segment. Variable productions, therefore, appear to be variable in their global shape, not just in their relationship to local targets.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript/Supplementary Files.

ETHICS STATEMENT
The study was reviewed and approved by the CUNY University Integrated IRB. Written and informed consent was obtained from the participant.

AUTHOR CONTRIBUTIONS
DW contributed conception and design of the study and wrote the first draft of the manuscript. W-RC organized the database, performed the statistical analysis, and wrote sections of the manuscript. DW and W-RC contributed to manuscript revision, read, and approved the submitted version.