Risk of Hearing Loss Caused by Multiple Acoustic Impulses in the Framework of Biovariability ()
1. Introduction
Hearing loss impacts about 360 million people worldwide (over 5 percent of the global population) [1]. One cause of hearing loss is exposure to loud sounds [2]. For example, impulse noise from explosive blasts and weapon firings can cause traumatic injuries to the auditory system [3] [4] [5].
In [4] hearing loss injury experiments were conducted at the University of California, San Diego, using chinchilla as the animal subjects. Each chinchilla was individually exposed to impulse noises generated by a small shock tube placed at various distances to the animal to cover a range of intensity levels; after the exposure, the injury status of the animal was measured and recorded; results from multiple individual animals were assembled into data sets to mimic the collective response of a crowd in responding to impulse noises. The main reason to choose chinchilla as a surrogate animal model for human exposures is that hearing capabilities are very similar between chinchilla and humans. However, the threshold for susceptibility of injury for chinchilla is found to be lower than humans. Thus, the chinchilla data were then shifted by 28 dBA, representing scaling from chinchilla to humans. Based on the scaled data, an empirical injury model was developed for human exposed to multiple sound impulses of equal intensity.
In our recent work [6] we use the framework of immunity to interpret the empirical dose-response relation from [3] for exposure to multiple sound impulses. More specifically, from the observed dose-response relation, we derived the (negative) synergistic effect of a sequence of impulses in causing injury. We revealed that the phenomenological effect of a preceding sound exposure event on the subsequent sound exposure event is always immunity. The study of [6] was based on the assumption that all subjects in the crowd are statistically indistinguishable and thereby it excluded the effect of biological variability. In this paper we focus on the effect of biological variability instead. This approach will provide an alternative view of the empirical dose-response relation from a different theoretical angle.
Biological variability (or “biovariability” in short) is defined as “the natural variability in a lab parameter due to physiologic differences among subjects and within the same subject over time” [7]. There are basically two types of biological variability, inter-individual and intra-individual. Inter-individual biovariability refers to “differences between subjects due to differences in diet, genetics or immune status” whereas intra-individual biovariability applies to “differences in the same subject over time due to diurnal cycles and other rhythms, biological repair mechanisms, dietary variables, aging, etc.” [7]. In this paper, we study the inter-individual biovariability. Prominent examples of inter-individual biological variability include diversities in biologically inherited traits such as height, weight and skin color. Due to biological variability there is always a range of responses displayed in measurements collected from a crowd of animal or human subjects. Therefore, it is reasonable to expect that some people are biologically less or more vulnerable to hearing loss injury than others. For example, one of the physiologic causes for variability is that body temperature is known to affect sensitivity to noise [3].
The presence of biovariability prompts us to wonder: what kind of role does the biological variability play when a crowd of subjects is theoretically exposed to multiple acoustic impulses? To address this question, we investigate the feasibility of explaining the observed logistic dose-response relation in the framework of biovariability. We mathematically derive the distribution density of individual injury probability in a crowd that will reproduce the observed dose-response relation. For several parameter values relevant for applications, we show that the derived distribution is a proper density function. Further, we study the asymptotic approximations of the density function and their usage in numerical computation of density function in finite precision arithmetic. Our mathematical results indicate that it is theoretically possible to interpret the observed logistic dose-response relation in terms of biovariability. This study on biovariability, together with the previous study on immunity effects [6] , offers two complementary views of the observed dose-response relation. These two very different views represent the two extreme theoretical boundaries for understanding the observed dose-response relation.
2. Mathematical Formulation
In experiments [3] , the injury risk of a crowd caused by a sound exposure was measured by conducting separate tests on animal subjects individually and assembling individual results to mimic the collective response of a crowd. A sound exposure event is a complicated process, described by the acoustic waveform, intensity, and duration. In [8] , several candidates were compared as an effective single metric for characterizing the overall injury causing effect of a sound exposure event. It was found that 8-hour Equivalent A-weighted Level (
) is the best predictor of injury [8] [9]. In the injury model developed in [3] , SELA (A-weighted Sound Exposure Level) was selected as the dose, which is mathematically equivalent to
, up to an additive constant. The observed injury risk follows the logistic dose-response relation:
(1)
For a crowd, the injury risk p is the average fraction of injured. In the logistic dose response relation (1), α is the shape parameter controlling the steepness of function, and
is the median injury dose. For a crowd, the median injury dose is the dose level at which 50% of the crowd is expected to be injured. We consider the hearing loss injury risk for a crowd of subjects with a wide spectrum of individual injury probabilities due to biovariability. That is, some subjects are biologically less or more susceptible to hearing loss injury than others. At the apparent median injury dose for such a crowd, 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 measured values of parameters α and
are respectively,
and
dB [3]. It was observed that parameter
remains the same for PTS injuries of all cut-off levels while the median injury dose
increases with the PTS cut-off level [3].
For a sequence of N identical impulses with SELA value = S, the effective combined dose is given by the dose combination rule [3] :
(2)
The injury risk of the crowd caused by the N impulses is
(3)
where parameters a and η are related to other parameters as
(4)
We express parameter a in terms of probability
by solving
,
(5)
This expression states that parameter a is the odds of injury of a hypothetical subject in the crowd with the average injury probability, which is the ratio of average injury fraction to average non-injury fraction of the crowd. In the subsequent discussion, we shall refer to parameter a simply as the odds of injury of an average subject. We shall avoid the term “average odds” since it is not the average of odds. Rather it is the odds for a subject with the average injury probability.
In this study, we assume that N impulses act independently from each other in causing injury, regardless of whether one is preceded by another (i.e., no immunity effect). Under the assumption of independence, we explore 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. When all events are independent of each other, the overall non-injury probability is simply the product of non-injury probabilities of individual events. Let
denote the intrinsic non-injury probability of a random subject in the crowd, in responding to one impulse. Here
is a random variable and we include symbol ω in notation explicitly to indicate the presence of biovariability. Let
be the probability density of random variable
. For one impulse, the average non-injury fraction of the crowd has the expression
For N impulses, the non-injury probability of an individual subject is
and the average non-injury fraction of the crowd is the average of
:
(6)
On the other hand, the injury risk follows the observed logistic dose-response relation (3), which relates to the average non-injury fraction as
(7)
Combining (6) and (7), we arrive at an equation for
, the distribution density of non-injury probability of a random subject in the crowd.
(8)
In the next section, we derive an analytical expression for the distribution density
based on Equation (8).
3. Analytical Results
We solve for density
in Equation (8). To proceed, we make a change of variables
. Since the non-injury probability
is constrained to interval
, variable
has the domain
. In terms of the new variable s, Equation (8) becomes
(9)
3.1. Power Series Solution
We expand the right-hand side (RHS) of Equation (9) as a power series of
.
(10)
Recall the mathematical identity
(11)
where the gamma function is defined as
We use identity (11) to map each power of
on the RHS of (10) back to a power of s.
(12)
Using this mapping, we express
on the left-hand side (LHS) of (10) as a power series.
(13)
Function
is the probability density of random variable
. The derivation above tells us that the theoretical prediction matches the observed logistic dose-response relation if and only if the density of
in the theoretical model has the power series in (13). It is not obvious at all whether or not function
defined by power series (13) is a proper density function. Below we first show that function
is well defined in its domain (
). That is, power series (13) converges for all values of
.
For large values of index k, the gamma function
is well approximated by the Stirling formula [10] [11] :
(14)
With the help of the Stirling formula, we calculate the radius of convergence using the ratio test:
(15)
It follows that function
is well defined by power series (13) for all values of s. To be a proper density function, however,
must satisfy two more properties:
(16)
and
(17)
We will rigorously demonstrate properties (16) and (17) for the cases of
and
. In addition, we will derive asymptotic approximations of
for large s, which are essential for practical computation of
in finite precision arithmetic. As we will see in the analysis, the proof approach differs quite significantly between these two cases. Thus, it is unlikely that the proof approach of either case can be directly extended to other values of η. Nevertheless, we conjecture that properties (16) and (17) are valid for all values of η.
To facilitate the analysis, we consider the cumulative distribution function (CDF) of random variable
:
. Integrating both sides of (13), we express
as
(18)
We use scaling
to get rid of parameter a. For simplicity and conciseness of presentation, we still denote the scaled variable by s and set
. After scaling, the cumulative distribution function
and the density function
have the expressions
(19)
(20)
Our general strategy is to derive a differential equation for
based on the power series around
. From the differential equation, we solve for
and
, and then derive their asymptotic approximations at
large s. Building on the theoretical insight gained in the cases of
and
, we then extend the asymptotic analysis to the case of
to derive an
essential numerical formula for practical computation of
in finite precision arithmetic.
3.2. The Case of
Setting
and differentiating (19) with respect to s, we have
(21)
which leads to a differential equation for
:
(22)
Using the integrating factor method and applying the condition
, we write the solution
as
(23)
It is straightforward to verify that
•
,
•
increases monotonically with s because
(24)
These two results lead directly to
and
.
Remark: It is fairly unusual for a solution of differential Equation (22) to converge to a constant level as s increases to
. In the homogeneous version of differential Equation (22), solution
exponentially diverges away from
as s increases. It is the unique combination of the condition at
and the RHS term in differential Equation (22) that makes
. If we change either or both of these two,
will diverge to infinity.
To derive the asymptotic approximations of
for large s, we change the integration variable to z defined by
. We write
as
(25)
For the integral in (25), at large s, the domain of dominant contribution is a small neighborhood near
. We expand
around
:
Each term in the power series of
gives us
Putting these results together, we write out the asymptotics of
for large s:
(26)
Differentiating with respect to s yields the asymptotics of
for large s:
(27)
Interestingly, after we concluded
, the asymptotics of
for large s can be derived in a simple way, directly from differential Equation (22), by assuming that
has an expansion of the form
(28)
Substituting expansion (28) iteratively into differential Equation (22) gives us
This set of equations yields the solution
The result of this simple method agrees with (26), which is rigorously derived without any additional assumption.
Remark: This simple method is quite powerful. By assuming the expansion form, we are able to derive the asymptotics of
at large s, based upon the power series around
.
Next, we study the case of
. We will see that the differential equation for
in the case of
behaves very differently from that in the case of
although the main conclusions regarding
and
remain the same.
3.3. The Case of
At
, differentiating (19) with respect to s leads to
(29)
which gives us a differential equation for
.
(30)
Using the integrating factor method and applying the condition
, we write the solution
as
(31)
In Appendix A, we prove that
given in (31) satisfies
1)
, and
2)
increases monotonically.
These two results imply that
is a mathematically proper density function.
To derive the asymptotics of
for large s, we use a change of variable
to rewrite solution (31) as
(32)
We expand
and
around
:
For large s, each term in the power series of
gives us
Substituting the power series into (32) and neglecting
terms, we write out the asymptotics of
for large s:
(33)
Differentiating with respect to s yields the asymptotics of
for large s:
(34)
Again, the asymptotics of
for large s can be derived directly from differential Equation (30) by assuming the expansion form (28). Below we use this approach to derive the asymptotics of
for large s in the case of
. We select the parameter value
to approximate the observed value
of
in experiments.
3.4. Asymptotics at Large s in the Case of
At
, differentiating (19) with respect to s leads to
(35)
which gives us a differential equation for
:
(36)
We derive the asymptotics of
for large s by assuming
•
satisfies
, and
•
has the expansion form (28) for large s.
Substituting expansion (28) iteratively into differential Equation (36), we obtain
(37)
Differentiating with respect to s yields the asymptotic of
for large s
(38)
3.5. Significance of Asymptotics at Large s for Computing
In this subsection, we discuss the necessity of using the asymptotic approximations in calculating density
for large s. The power series of
around
, given in (20), has the nice property that it converges analytically for any s. In practical computation with finite precision arithmetic, however, the numerical convergence of (20) is limited by the finite precision and the limitation is fatal at large s. To see this limitation, we examine the net sum of the power series and compare the net sum with individual terms in the power series. When the largest term is 1016 times the net sum or more, implementing power series (20) with double precision arithmetic will lose all numerical accuracy in computing the net sum.
To simplify the discussion, we write power series (20) in a concise and abstract form
(39)
We study the case of
as an example. The power series converges for
any s and the exact net sum is
, which has the asymptotic approximation (38) at large s. We use the leading term in (38) to estimate the net sum for large s:
(40)
This gives us the general magnitude and trend of
. In particular,
decreases to zero as s increases. In power series (39), we write the absolute value of the k-th term as
(41)
The largest term in power series (39) is found by maximizing function
with respect to index k. For mathematical convenience, we approximate this discrete maximization as a continuous maximization process. We maximize
with respect to z as a continuous variable. Using the Stirling approximation (14), we write
as
(42)
We maximize
with respect to z as a continuous variable. The derivative of
is
(43)
The derivative indicates that for large s, the maximum of
is attained
approximately at
. Consequently,
is a fairly tight lower
bound on
. That is, the true maximum may be a bit larger but not
too much larger than
:
(44)
The largest term (in absolute value) in power series (39) is approximately
(45)
In summary, for large s, we estimated the net sum and the largest term in power series (39):
• Net sum:
,
• Largest term:
.
The ratio of largest to net sum is
(46)
At
, for example, the largest term in (39) is approximately
; the net sum of (39) is
; and the ratio of largest to net sum is
. Recall that the machine precision of the IEEE double precision is about
. Clearly, at
, computing function
by numerically summing terms in power series (39) will not yield any accuracy when implemented in double precision arithmetic. In addition, as s increases, the ratio of largest to net sum, given by Equation (46), grows exponentially with s. As a result, adopting a higher precision arithmetic will only enable us to accommodate a slightly larger range of s. For example, at
, the ratio of largest to net sum is
which will completely wipe out the numerical accuracy of quadruple precision arithmetic. In conclusion, although power series (39) converges theoretically for any s, it is not a practical numerical method for computing function
at large s. For large values of s, the asymptotic approximation (38) is absolutely essential in providing a viable way for computing function
.
Let us calculate
, the index value at which the largest term occurs in power series (39). For large s, index
satisfies approximately
(47)
Solving for
from this equation, we obtain
For
and
, we have
.
Figure 1 plots the k-th term of power series (39) vs index k for
and
. The largest term occurs at indices k = 282 and k = 283, which correspond very well to the predicted value of
.
Figure 1. Plot of
vs k for
and
. The inset shows details of the plot near the location of maximum.
It is worthwhile to find
, the cut-off index beyond which all terms in the power series are smaller than a threshold value specified by a relative error tolerance (ε). Specifically, we find index
such that
(48)
where ε is the relative error tolerance. With the expression of
given in (41), we cast it approximately into a continuous problem
Applying the Stirling approximation (42) on the LHS and applying the leading term asymptotics of
for large s (40) on the RHS, the constraint on
becomes
(49)
For large s, the largest term of power series (39) occurs at
. We expect the cut-off threshold
to be a small multiple of
, which suggests us to re-write the inequality for
approximately in terms of
:
(50)
Function
is well approximated by its quadratic expansion.
Using this quadratic approximation, we solve for
from inequality (50).
(51)
The corresponding index
is
(52)
For
,
, and
, we have
.
Note that summing the first
terms in the power series (39) is only a necessary condition for achieving a relative accuracy of
in calculating
. It is not a sufficient condition. In the numerical computation, the result is polluted by the accumulation of round-off errors caused by finite precision arithmetic. Thus, the numerical accuracy is limited by the machine precision of the number representation system multiplied by the largest term in the power series. As we showed above, with double precision arithmetic, even at a moderate
, the round-off errors completely wipe out the numerical accuracy.
3.6. Density of Individual Injury Probability in the Crowd
In this subsection, we write out
, the distribution density of individual injury probability, which shows the biological variability of the crowd in its susceptibility to hearing loss injury.
In all analysis above, function
is the density of
after scaling
. Density of
before the scaling is given by
(53)
Going back to non-injury probability
, we write out the density of individual non-injury probability
in terms of
.
(54)
The density of individual injury probability
is
(55)
where
is the density of the post-scaling random variable
. Function
is calculated as follows.
• For
,
is numerically computed using the partial sum of first
terms in the power series around
given in (20).
(56)
Here
is the number of terms to include in summation to make the truncation error smaller than 10−16.
is estimated using (51) and (52). Since the estimate of
is not very reliable for small s, we use at least 500 terms in summation.
• For
,
is evaluated using the first
groups of terms in asymptotic approximation (38) for large s. In (38), only the first 3 groups are explicitly shown. Now let us write out a general formula for the partial sum of the first
groups of terms in the asymptotic approximation. The general formula is in the form of double summation and is suitable for numerical implementation.
(57)
The threshold (
) for switching from power series around
to asymptotic approximation at large s is discussed in next section.
We need to point out that in asymptotic approximations (38) and (57), at a fixed value of s, if we include more and more groups of terms in the approximation, eventually the sum diverges. This is caused by the divergence of infinite series
In this infinite series, the gamma function coefficient increases much faster than exponential growth. At any fixed value of s, as index n increases, the n-th
term behaves like
, which diverges to infinity. At a fixed and moderate
value of
, however, asymptotic approximation (57) provides a vital numerical tool for computing function
at large s when the power series around
loses numerical accuracy due to accumulation of round-off errors.
4. Numerical Results and Discussion
In this section, we study the numerical behaviors of the power series around
and the asymptotic approximations at large s. Based on the numerical results, we develop a numerical procedure for computing
that is accurate for all values of s, small or large. The procedure is already mentioned in the previous section and is based on switching between power series (56) and asymptotics (57). Using this numerical procedure, we examine the distribution density
of the individual injury probability in the crowd. We also plot the cumulative distribution function (CDF) of the injury probability because the CDF illustrates the effect of parameter a on the distribution of individual injury probability. Here parameter a is the observed odds of injury of an average subject in the crowd.
All simulations in this section are for the case of
, which is meant to
approximate the observed value of
in experiments. We first study the errors in numerical implementation of power series (56) in IEEE double precision. When sufficiently large number of terms are included in the power series summation to keep the theoretical truncation error well below the round-off errors, the numerical errors are mainly caused by accumulation of round-off errors. To calculate the numerical errors, we need a very accurate solution. For that purpose, we implemented all operations and functions (including
,
and
) in a very high numerical precision by representing each number using an array of 200 integers to achieve 1600 significant decimal digits. Numerical error is calculated as the difference between the numerical solution in double precision and the high precision solution. Figure 2 plots the absolute error and relative error vs s when the power series around
is implemented in IEEE double precision. In the previous section, we derived that the largest term in power series (56) is
. Based on that, we expect the absolute numerical error to grow
approximately proportional to
. In Figure 2, the gray dashed line is a fitting of the form
to the numerical absolute errors, confirming the theoretical prediction.
The exponential growth of the accumulated round-off errors with s implies that, for large s, the numerical summation of power series (56) is not a good tool for computing
. For large s, we turn our attention to asymptotics (57) of
. The numerical errors in using asymptotics (57) to calculate
are mainly due to the approximation errors at finite s. In the previous section, we showed that at any fixed value of s, if we include more and more terms in the asymptotics, it eventually diverges. At a finite s, the approximation error of asymptotics cannot be made arbitrarily small by including more terms even if we work in exact arithmetic (i.e., no round-off error). The approximation errors can only be reduced by increasing s. Thus, asymptotics (57) provides a good approximation of
only when we use a fixed number of terms and when s is at least moderately large. Figure 3 displays the approximation errors in asymptotics (57) vs s when respectively 6 groups of terms and 12 groups of terms are used in approximation. For
, the 12-group approximation is more accurate while for
, the 6-group approximation is more accurate. With a fixed number of groups of terms, the approximation is more accurate as s is increased.
Figure 2. Numerical errors vs s when power series (56) of
is implemented in IEEE double precision. Here the numerical errors are caused by the accumulation of round-off errors. The errors grow exponentially with s because the largest term in summation grows exponentially. The gray dashed line is a fitting of the form
.
These numerical findings agree with the theoretical predictions in the previous section.
From the results of Figure 2 and Figure 3, we see that a good numerical strategy for calculating
is to use the power series around
when s is not too large and switch to the asymptotics when s is large enough. Figure 4 compares the errors of these two formulas: 1) the numerical error when implementing power series summation (56) in IEEE double precision, and 2) the approximation error in asymptotics (57). The accumulation of round-off errors in power series grows with s while the approximation error in asymptotics decays as s is increased. Figure 4 suggests
as an optimal threshold for switching from power series (56) to asymptotics (57). When this switching strategy is implemented in IEEE double precision, it provides a numerical procedure for calculating
that is accurate to 10−8 for all values of s in
.
With the switching mechanism developed above for calculating
, we study the distribution density and cumulative distribution function of the individual injury probability. In Figure 5, we plot density
of the individual injury probability p in the crowd. Density
is given in equation (55), expressed in terms of function
. Figure 5 indicates that density
diverges to infinity at
and
. Approximately, distribution
can be viewed as the mix of 3 component distributions:
1) a δ-function distribution at
,
2) a nearly uniform distribution in
, and
3) a δ-function distribution at
.
Because of the divergence of density at
and at
, the weights of the two δ-function components cannot be readily read out from the plot of density
Figure 3. Approximation errors in asymptotics (57) vs s when respectively 6 groups of terms and 12 groups of terms are used in approximation. With each fixed number of groups of terms, the asymptotic approximation becomes more accurate at larger s.
Figure 4. Comparison of two errors: 1) the error in numerical summation of power series around
, and 2) the approximation error in asymptotics for large s. The intersection of the two suggests a threshold for switching from one to the other.
. To see the weights, we examine
, the cumulative distribution function (CDF) of p, which clearly illustrates the weight of a δ-function component as a jump in function value of the same magnitude.
The CDF,
, has the expression
(58)
where
is CDF of the post-scaling random variable
, expressed as a power series around
in (19). While theoretically, power series (19) converges for any s, the practical numerical convergence is limited
Figure 5. Density
of the individual injury probability p for
where parameter a is the observed odds of injury of an average subject in the crowd. The logarithmic scale is used in the vertical direction because
diverges to infinity at
and
.
by the precision of the number representation system. In IEEE double precision, we adopt a switching numerical procedure for calculating
, similar to that for calculating
.
• For
,
is numerically computed using the partial sum of first
terms in the power series around
given in (19).
(59)
• For
,
is evaluated using the first
groups of terms in the asymptotics for large s given in (37). We rewrite it as a general formula below.
(60)
Figure 6 compares plots of
respectively for
,
, and
. Parameter a is the observed odds of injury of an average subject in the crowd. As parameter a increases from small to one to large, the δ-function component at
decreases and the δ-function component at
increases by the same amount while the distribution in between is approximately unchanged.
5. Conclusion
In this study, we examined the observed logistic dose-response relation for hearing loss injury caused by exposure to multiple acoustic impulses, in the mathematical framework of biovariability. Previously, we interpreted the
Figure 6. Plots of
, CDF of the individual injury probability, respectively for
,
, and
. Parameter a is the observed odds of injury of an average subject in the crowd. As parameter a increases, its most prominent effect is to shift F(p) down nearly uniformly while keeping
and
.
observed dose-response relation based solely on the immunity effect under the assumption that all individual subjects are statistically identical (i.e., no biovariability). In the current study, we view the problem from a completely different angle; we explored the possibility of understanding the observed dose-response relation based solely on the biovariability of the crowd under the assumption that all sound exposure events act independently from each other in causing injury regardless of whether one event is preceded by another. We found that theoretically the biovariability alone is indeed capable of explaining the observed dose-response relation. We derived an analytical expression for the distribution density of individual injury probability in the crowd that produces the observed results. We also derived asymptotic approximations to the distribution density and the cumulative distribution function. The asymptotic approximations make it possible to compute the distribution density in finite precision arithmetic, in parameter regions where the analytical expression suffers catastrophically from loss of accuracy due to the accumulation of round-off errors. A robust numerical procedure for calculating the distribution density was developed based on switching between the analytical expression and the asymptotics. The resulting numerical procedure is accurate in all parameter regions, providing a practical tool for computing the distribution density. The mathematical framework constructed in the current study, along with the theoretical results and the numerical tools obtained, paves a pathway for further investigating the presence and effects of biovariability.
6. Disclaimer
The authors would like to 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.
Appendix A
We prove that
given in (31) satisfies
1)
, and
2)
increases monotonically.
Item 1) is obtained directly by applying L’Hospital’s rule.
(61)
To prove 2), we introduce function
We write
in terms of
and define function
:
(62)
We only need to show
. The derivative of
has the expression
(63)
From the definition,
and
satisfy
(64)
(65)
where
The sign of
tells us that
(66)
In other words,
attains its minimum at
. Thus, to show
, we only need to show
. We use integration by parts to rewrite
as
(67)
Function
satisfies
Using the sign of
, we have
(68)
Combining these results, we conclude that
and thus, function
increases monotonically.