Sound is an indispensable part of our life and we experience sound every day. A common way to measure the amount of sound is the decibel (dB)  . Sounds of less than 75 dB are at safe levels that do not damage our hearing. However, any sound above 85 dB is potentially harmful and can cause hearing loss. Examples of harmless sounds are normal conversation (60 dB), the humming of a refrigerator (45 dB) whereas harmful sounds include noise from lawn mowers (90 dB) and gun shots or firecrackers (both 150 dB)  . The risk of hearing damage depends on the power of the sound as well as the length of exposure.
Hearing loss is a common health problem among veterans. In order to protect warfighters, starting 1960s the US Army conducted and funded research to assess the risk of hearing loss caused by intense impulse noise from explosive blasts and weapon firings  . Recently Dr. Chan and his collaborators  developed a dose-response model for the assessment of injury caused by impulse noise and a model for the possible recovery afterwards, based on chinchilla data. Chinchillas share similar hearing capabilities as humans and thereby are commonly used for hearing-related experiments.
In  , we interpreted the empirical dose-response relation from  for exposure to multiple sound impulses in the framework of immunity. In  , we viewed the empirical dose-response relation from a completely different angle, in the framework of biovariability. Together in these two studies   , we showed that it is possible to interpret the empirical dose-response relation from either of the two extreme cases: immunity or biovariability. Here we would like to further our study in  to demonstrate that the derived distribution density of injury susceptibility in  is well-posed.
2. Mathematical Formulation of the Problem
In experiments  , the injury risk of a crowd caused by a sound exposure event is described by the logistic dose-response relation:
Here the dose of the sound exposure event is defined by the SELA (A-Weighted Sound Exposure Level) in units of dBA. The injury risk p of the crowd represents the average injury fraction of the crowd. In the logistic dose response model (1), the parameter α determines the steepness of function p while D50 denotes the median injury dose. For a crowd, the median injury dose is the dose level at which half of the population is expected to be injured. For a crowd of subjects with a wide spectrum of heterogeneous individual injury probabilities, at the apparent median injury dose, a particular subject's individual injury probability may be below or above 50% due to biovariability. For injury of permanent threshold shift (PTS) > 25 dB, the values of parameters α and D50 are found to be and dB, respectively  . It was also noticed that the parameter value α remains unchanged ( ) for PTS injuries of all cut-off levels whereas the median injury dose D50 rises with the PTS cut-off level  .
In the framework of the logistic dose-response relation, the injury risk of the crowd caused by a sequence of N acoustic impulses is given by the expression
where is the effective combined dose for the whole sequence of impulses as a single sound exposure event. For a sequence of N identical impulses each with SELA value S, the effective combined dose, , was observed to follow the dose combination rule  :
Thus, for a sequence of N impulses each with SELA value S, the injury risk takes the form
where parameters a and η are defined as
Parameter a is related to probability by
That is, parameter a is the injury odds of a hypothetical subject in the crowd with the average injury probability, , in responding to a single acoustic impulse.
In  under the assumption that N acoustic impulses act independently from each other in causing injury, regardless of whether one is preceded by another (i.e., no immunity effect), we explored the possibility of interpreting the observed logistic dose-response relation for a crowd in the framework of biovariability. For mathematical convenience, we consider non-injury probability instead of injury probability. Let denote the individual non-injury probability of a random subject in the crowd, in responding to one acoustic impulse. Here is a random variable, due to the presence of biovariability. Let be the distribution density of random variable . Mathematically in the framework of biovariability, the average non-injury fraction for N acoustic impulses is expressed in terms of as
On the other hand, experimentally, the injury risk was observed to follow the logistic dose-response relation (4), which relates to the average non-injury fraction as
For the theoretical model of biovariability to reproduce the experimentally observed results, the distribution density has to satisfy an equation obtained by combining (7) and (8), which we write out below:
In  we solved Equation (9) analytically by constructing a power series
expansion in new variable . Since the non-injury probability is constrained to interval , variable has the domain . The distribution density of random variable is
For density , Equation (9) becomes
In  we derived a power series solution for :
In power series (12), the gamma function in the denominator of the coefficient grows faster than any exponential function of k. As a result, power series (12) converges for all values of s. It follows that function is well defined by power series (12) for all values of s. To be a proper density function, however, must satisfy the two properties below:
In  we rigorously proved these two properties for the special case of and the special case of . The analysis procedure differs quite significantly between these two special cases. It is highly unlikely that the particular analysis approach used in either of these two special cases can be directly extended to the general case of arbitrary η. In the current study, we aim at verifying semi-analytically the two properties for any given value of η. For that purpose, we need the numerical capability of calculating function for all values of s, from small to large. Power series (12) has the nice property that theoretically it converges for all values of s. In practical computations, however, at large values of s, the numerical accuracy of the power series summation is completely wiped out by the accumulation of round-off errors. In the power series summation, as s increases the net sum decreases while the magnitude of the largest term grows exponentially with s  . The combined effect of these two factors magnifies catastrophically, at large s, the influence of round-off errors on the numerical accuracy of the net sum. For practical computation of in finite precision arithmetic, we need a robust numerical formula of at large s. In the next section, we derive a general asymptotic approximation of at large s. The synthesis of the power series and the asymptotics will provide a practical numerical tool for computing the distribution density for all values of s at any given value of parameter η.
3. Asymptotics of g(s) at Large s
Now we derive a general asymptotic approximation of at large s when η is a rational number. We then reasonably conjecture that the same asymptotic approximation is also valid even when η is irrational. In practical computation of , the case of irrational η actually does not apply since all numerical calculations are carried out in finite precision arithmetic, using only rational numbers.
A rational number η takes the form where both m and n are positive
integers. We rewrite the power series of in terms of the reciprocal gamma function as follows:
where is the reciprocal gamma function defined as
The advantage of working with the reciprocal gamma function is that it is well defined and is analytic everywhere. In comparison, the gamma function diverges at all non-positive integer values of z. The reciprocal gamma function has the property
Using this convenient property when differentiating , we have
Differentiating repeatedly m times, we obtain a differential equation for .
We derive an asymptotic approximation of at large s based on this differential equation. We proceed with the assumption that converges to 0 as s goes to . This assumption is reasonable although not directly derivable from the power series form of . Under this assumption, the asymptotics of at large s takes the form
Differentiating the asymptotic form m times, we get
In (21) the first term is of the order while all other terms are asymptotically smaller than . Using the result of (21), we rewrite (19) as
which leads to
Equating the leading terms on both sides of (23) yields . Substituting this value of back into (23) gives us
Notice that in (24), the smallest term in the summation is which occurs at . Thus, all terms in the summation are indeed asymptotically larger than .
In summary, expression (24) gives an asymptotic approximation of at large s, accurate up to the order . Recall that (24) is derived by differentiating each of power series (15) and asymptotic form (20) m times, and then combining the results. To derive a general asymptotic approximation, we differentiate each of the power series and the asymptotic form times where r is a positive integer. In the above, asymptotics (24) is the outcome in the special case of . In the general case, differentiating power series (15) times yields.
Differentiating asymptotic form (20) times, we get
where we have used . Combining (25) and (26), we obtain
Since integer r can be made as large as we like, in (27) we simply use the infinite series as a symbolic general asymptotics, with the understanding that a particular asymptotic approximation will use only a partial sum of the infinite series. We need to point out that the infinite series in (27) actually diverges for all values of s. So the infinite summation does not have a well defined sum. Instead, the infinite series serves only as a symbolic general asymptotics. Summation of moderate number of terms, however, will provide an accurate asymptotic approximation of at moderately large s and beyond. We will examine the approximation errors in details later. In the analysis above, we did not assume that integers m and n are prime to each other. It turns out that asymptotics (27) can be written in terms of η only, without any reference to m or n. The expression in terms of η gives us the general asymptotics, which does not depend on a particular rational form of η:
The asymptotic expansion (28) depends only on η and is invariant with respect to different rational representations of η. It is plausible to conjecture that asymptotics (28) is also valid for irrational η although it is derived in the case of rational η. This conjecture cannot be tested numerically since all computations use a finite precision number representation system, which is a subset of all rational numbers. In the next section, in preparation for the numerical verification of properties (13) and (14), we build the necessary numerical tools for computing function .
4. Accurate Evaluation of g(s) in Finite Precision
In this section, we develop a practical numerical method for computing in IEEE double precision, over the whole range of s and at any given value of parameter η.
Function is defined straightforwardly by power series (12). Theoretically can be evaluated as accurately as we like by including sufficiently large number of terms in summation and carrying out the computation in arithmetic of sufficiently high numerical precision. Practically, with the IEEE double precision arithmetic, the numerical accuracy of the power series summation, at large s, is completely ruined by the round-off errors from terms of the largest magnitude. In  we showed that the largest term grows roughly exponentially with s and it has the behavior.
Even at a moderate large value of , the largest term in the summation is more than 1016, which in general will pollute the numerical value of with an error of magnitude 1 or bigger. Thus, at large s, the power series summation is not a workable numerical tool for accurately calculating in finite precision arithmetic.
The infinite series in the asymptotic expansion (28) diverges for all values of s. As a result, it does not make sense to include in the asymptotic approximation a very large number of terms from (28). When a moderate number of terms are used, however, the partial sum of (28) provides an accurate approximation of at moderately large s and beyond. For a fixed number of terms, the larger the value of s is, the better the approximation. Therefore, at large s, function can be evaluated fairly accurately by employing an asymptotic approximation with a suitable number of terms.
The contrast and complementary behaviors of the power series around and the asymptotics at large s suggest that a viable numerical strategy is to use the power series summation for small s and switch to the asymptotic approximation when s is above a threshold ssw, which is yet to be specified. The success of this numerical strategy depends on that there is an overlapping region of intermediate s in which both the power series summation and the asymptotic approximation will yield reasonably good accuracy. Without this intermediate region, if the valid region of the power series summation is separated by a gap from the valid region of the asymptotic approximation, cannot be calculated accurately in the gap region. The existence of an overlapping region of intermediate s also provides us a numerical mechanism for identifying the overlapping region and selecting an optimal threshold ssw for switching from one numerical formula to the other.
In the overlapping region, both the power series summation and the asymptotic approximation are reasonably accurate. Accordingly, the difference between the two numerical formulas should be fairly small inside the overlapping region. Below or above the overlapping region, only one of the numerical formula is very accurate while the other is not. Consequently, outside the overlapping region, the difference between the two will be significantly more pronounced than inside the overlapping region. To identify the overlapping region, we examine the difference between the power series summation and the asymptotic approximation as s increases from small values to large values. The magnitude of the minimum difference indicates the existence (or non-existence) of the overlapping region; the location of the minimum difference suggests an optimal threshold ssw for switching. To proceed along this line, we introduce two notations.
• = power series summation (12)
• = asymptotic approximation (28) with terms up to ...
Throughout this paper, all numerical results are computed in IEEE double precision arithmetic. In this section, simulations are focused on the case of , the observed value of parameter η in experiments. We will explore other values of η in subsequent sections.
In Figure 1, we plot the difference between and as a function of s for several different values of Ng. Three asymptotic approximations respectively with , and are tested. For all 3 values of Ng, especially for and , Figure 1 demonstrates clearly the existence of an overlapping valid region for the two numerical formulas. For s values smaller than 17, there is a visible discrepancy among the three curves because in this region is primarily attributed to the approximation error in which depends heavily on Ng when s is not very large. For s values bigger than 17, the three curves almost coincide with each other due to the dominant effect of round-off errors in which is independent of Ng. The overlapping region is in a neighborhood around   . The asymptotic approximation with has the best performance since it reaches the lowest minimum difference and attains the minimum difference at a smaller value of s, indicating that it is already valid even when s is not very large. In our subsequent simulations, we shall select Ng using this strategy. For , Figure 1 suggests that an optimal threshold ssw for switching from power series summation to asymptotic approximation is about . Based computing function .
Figure 1. Difference between the power series summation and the asymptotic expansion for various values of Ng.
• For , is computed using the partial sum of the first Kf terms in power series summation (12).
on these numerical findings we adopt the numerical procedure below for where Kf is the number of terms needed to make the truncation error well below the machine precision of IEEE double precision. In computations, we make the truncation error smaller than 10−20. Theoretically, Kf has an a priori estimate expressed in terms of s when s is moderately large. In practice, Kf is determined automatically in the numerical summation process by monitoring the magnitude of terms.
• For , is calculated using the partial sum of terms up to in the asymptotic approximation (28).
The choice of and above is based on numerical minimization of with respect to in the case of . This particular set of is for computing function at . In a similar fashion, at each different value of η, an individual set of is determined for that η and then used in evaluating .
In the next section, we apply the numerical procedure described above to verify properties (13) and (14) numerically, and thus demonstrate the well-posedness of the distribution density.
5. Numerical Verification of Well-Posedness
We first verify that defined by power series (12) is positive for all values of at any parameter value of . To examine both the sign and the magnitude of , we use mapping , defined below, to display as . Let
where is a design parameter depending on the range we like to focus on.
The mapping has several design features for showing the sign and for accommodating a huge range over many orders of magnitude.
• is an odd function of z, clearly showing the sign of z.
• is a monotonically increasing function of z, preserving any trend of z.
• When is significantly below , the mapping displays z in a linear scale:
• When is significantly above , the mapping displays z in a logarithmic scale:
We calculate vs s numerically for several representative values of , and plots in Figure 2 with . Function is positive for all values of η we examined.
Next we verify that for parameter . We integrate the power series (12) to write out the cumulative distribution function (CDF), which becomes
Again, theoretically power series (33) converges for all s, making a well defined function for all s. But in numerical computations with IEEE double precision, at large s, power series summation (33) suffers catastrophically from complete loss of accuracy. As a result, using the power series summation to compute at large s is not viable for demonstrating . Instead, we consider the complementary cumulative distribution function at large s, defined as
Figure 2. Plots of for several values of parameter η. The mapping defined in (32) is designed for showing the sign and for accommodating a wide range of quantity z. The plots demonstrate that function is positive for all values of parameter η tested.
To verify , we only need to demonstrate numerically that at some value of s. This allows us to select a value of s such that both and can be computed accurately.
For , can be accurately calculated using a partial sum of power series (33)
For , is well approximated by asymptotics (28). Using a partial sum of (28) with terms up to to replace in the integral of , we write out an asymptotic approximation for
Results (35) and (36) suggest that quantity has the optimal numerical accuracy if we evaluate it at . Thus, we compute quantity and use it to judge if is satisfied.
Figure 3 plots vs parameter η. It is clear that for any parameter value η in , the assertion is indeed valid within the numerical approximation error.
Figure 3. The difference between and 1 for different values of η.
In summary, we have numerically verified 1) for and 2) . It follows that function defined by power series (12) is mathematically a proper distribution density.
With the property established, we write out a unified numerical procedure for computing function over the full range of ,
For readers’ convenience, we also summarize below the unified numerical procedure for computing function over the full range of ,
6. Pade Approximations
In the previous section, we verified . We now impose this property as a constraint at and construct a Pade approximation   for based on its power series around . As we will see, the Pade approximation provides an accurate and efficient approximation over the full range of .
Note that one can set and start the summation at in (39). Taking these features of function into consideration, we adopt a Pade approximation of order , of the form
where follows from , and follows from and normalization. There are unknown coefficients in Pade approximation (40). To determine these coefficients, we multiply both (40) and (39) by , and then match the terms for . The product of two power series has the expression:
For k in the range of , equating the corresponding coefficients of terms on the left-hand side and on the right-hand sides yields equations for the unknowns .
where is known and has the expression
Equation (41) is an linear system for . Coefficients are determined by solving linear system (41). Once coefficients are known, we write out coefficients by matching the coefficients of terms for .
To estimate the error of Pade approximation defined in (40), we use computed with the unified numerical procedure (37) as the “exact” solution to compare with. We calculate the difference between the numerical value of and the Pade approximation . Figure 4 shows vs s for parameter value . Four Pade approximations, respectively with n = 3, 4, 5, and 6, are shown where n is the highest power used in Pade approximation (40). For n = 4, the approximation error of is already below 10−8, which is similar to the errors of both the power series summation and
Figure 4. Discrepancy between and Pade approximation .
the asymptotic approximation in the overlapping region around . The numerical error of procedure (37) varies with the magnitude of s. The numerical error is the largest in the overlapping region: below the overlapping region, the power series summation is less polluted by the round-off errors and thus is more accurate; above the overlapping region, the asymptotic approximation becomes more accurate. The numerical accuracy of calculated using (37) is significantly higher than 10−8 when s is outside the overlapping region. This property of will help us decipher the error behavior in higher order Pade approximations.
When n is increased to , the difference is below 10−10 outside the overlapping region, implying that the error of Pade approximation is also below 10−10 outside the overlapping region. The difference increases significantly in the overlapping region. However, it is highly unlikely that the approximation error of Pade approximation jumps significantly only in the overlapping region while remaining below 10−10 outside the overlapping region. The Pade approximation consists of one rational function for all values of s; it does not involve any switching. It is much more likely that the approximation error of Pade approximation actually remains below 10−10 over the full range of s; the significant increase in is solely caused by the increased numerical error of in the overlapping region. If this is true, then for , the Pade approximation is already more accurate than the unified numerical procedure (37) in IEEE double precision. The smaller numerical error of the Pade approximation is mainly attributed to that it has only a few terms, and subsequently, is much less affected by round-off errors in IEEE double precision. For , Figure 4 shows that the increase of near the overlapping region is much more pronounced than in the case of . The pattern of increase strongly suggests that it is caused by the increased numerical error of near the overlapping region. Figure 4 indicates that the true numerical error in Pade approximation is very likely below 10−12 throughout the full range of s, much more accurate than the unified numerical procedure (37) in IEEE double precision.
We carry out single precision computations to support the assertion we made above that in finite precision arithmetic, Pade approximations can be more accurate than the power series even though the power series is theoretically exact in infinite precision arithmetic. We use the double precision result of as the exact solution to compare with. We compute power series summation, asymptotics, and Pade approximations in single precision, and then examine the numerical errors in single precision results. Figure 5 shows the error behaviors of single precision results. As s increases, the power series summation starts losing accuracy due to the exponential growth of the largest term and the associated round-off error in summation. Meanwhile, as s increases, the approximation error in asymptotics decreases and its numerical accuracy improves. In contrast, the numerical errors in Pade approximations remain fairly steady with respect to s and decays very rapidly as n is increased. Pade approximation is already significantly more accurate than both the power series summation and the asymptotics in a large neighborhood of the overlapping region (for single precision arithmetic, the overlapping region is around ). It is evident in Figure 5 that the numerical error of is primarily caused by round-off errors and its true approximation error is below the machine epsilon of single precision
Figure 5. Errors of power series, asymptotics and Pade approximations in single precision computations.
system. In single precision, the actually realized numerical accuracy of is uniformly much higher than that of both the power series summation and the asymptotics over the full range of s.
Next, we construct a Pade approximation for density . Power series of around has the form
Note that since , we have in
(44), which suggests us to adopt a Pade approximation of order , of the form
There are 2n unknown coefficients in the Pade approximation of . To determine the coefficients, we multiply both (45) and (44) by , match the coefficients of terms for to form a linear system, and then solve for the unknowns, following a procedure similar to the one used in the Pade approximation of .
To assess the error of Pade approximation given in (45), we use computed with the unified numerical procedure (38) as the “exact” solution to compare with. Figure 6 shows vs s for parameter value . Four Pade approximations, respectively with, n = 3, 4, 5, and 6, are shown where n is the highest power used in Pade approximation (45). The behaviors of the Pade approximations for density are similar to those for the CDF . In IEEE double precision, the actually realized numerical accuracy of the Pade approximations with and is significantly better than that of numerical procedure (38). Again, in IEEE double precision, the smaller numerical error of the Pade approximations is mainly attributed to the fact that it contains only a few terms, and its numerical results are much less contaminated with round-off errors.
We studied the biovariability of a crowd for hearing loss injury, in the form of heterogeneous injury susceptibility. We constructed a unified numerical procedure for computing the distribution density of injury susceptibility that reproduces the observed logistic dose-response relation in a crowd. The unified procedure combines the advantage of power series expansion for small values of argument and the advantage of asymptotic approximation for large values of argument. It switches between these two approaches to achieve a numerical
Figure 6. Differences between density and its Pade approximations .
accuracy of 10−8 or better with IEEE double precision, over the full range of argument. Using this unified procedure, we verified numerically that for all parameter values, the derived distribution density, i) is non-negative everywhere and ii) integrates to one. These results establish numerically that the derived distribution is indeed a proper density for all values of parameter, and thus, is well-posed. Furthermore, we developed efficient and accurate Pade approximations for the distribution density and for the cumulative distribution function. In the computational environment of IEEE double precision, Pade approximations actually yield a much higher realized numerical accuracy than that of both the asymptotic approximation for large argument value and the power series for small argument value. The superior performance of Pade approximations is attributed to the fact that it attains high theoretical accuracy with only a few terms, which leads to less contamination with round-off errors and better realized numerical accuracy. In conclusion, we verified numerically that the observed logistic dose-response relation can be explained solely based on a valid distribution of injury susceptibility. Rigorous proof of the well-posedness of the derived distribution density, however, still remains open.
The authors thank the Joint Non-Lethal Weapons Directorate of US Department of Defense for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the US Government.
 Murphy, W.J., Khan, A. and Shaw, P.B. (2011) Analysis of Chinchilla Temporary and Permanent Threshold Shifts Following Impulsive Noise Exposure.
 Wang, H., Burgei, W.A. and Zhou, H. (2017) Interpreting Dose-Response Relation for Exposure to Multiple Sound Impulses in the Framework of Immunity. Health, 9, 1817-1842.
 Wang, H., Burgei, W.A. and Zhou, H. (2018) Risk of Hearing Loss Caused by Multiple Acoustic Impulses in the Framework of Biovariability. Health, 10, Article ID: 84786.