The problem of finding the distribution of self-normalized sum of a Cauchy-generated sample has been considered since at least 1969 (see, for example,  ), but solved (by proposing a concrete numerical algorithm) only in 1973 by  , our key reference. After that, the research has been focused mainly on proving general statements concerning self-normalized sums, without dealing with specific distributions such as Cauchy—see, for example,  .
In this article, we demonstrate how characteristic functions (CHF) are used in Statistics to find a distribution of a specific sample statistic (a function of individual observations). We do this by using a single comprehensive example, namely finding the limit of the probability density function (PDF) of
where is a random independent sample (RIS) of size n from a Cauchy distribution with zero median. This goal has already been achieved (in a more general setting) by  ; our main (rather pedagogical) reason for extending their presentation is to make it accessible to graduate and advance undergraduate students.
First we recall that the CHF of a random variable X is defined by
where stands for the X’s PDF, for the corresponding expected value, and i is the purely imaginary unit; note that the real part of is always an even (an alternate name is symmetric) function of s, while the purely imaginary part is odd. Also note that when is even, the corresponding CHF must be real (its purely imaginary part becomes zero because the integrand is odd).
(where is the joint PDF of X and Y) defines the joint CHF of two (not necessarily independent) random variables X and Y; when they are independent, we have
(product of the individual CHFs). This implies that (when independent), the CHF of is given by , further implying that the CHF of a sample mean of n independent is given by
When two random variables are combined into one (e.g. defining ), the CHF of U is given by
Under very general conditions, knowing a generating function enables us to reverse the process and find the corresponding PDF thus:
To learn more about the mathematical details of transforming PDF into CHF and back the reader may like to consult  .
2. Joint CHF and Its Limit
Assume a RIS of size n from a Cauchy distribution with the following PDF
The joint CHF of X and (seen as distinct random variables) is then given by
This implies that the joint CHF of
is (replace s by and t by in (9), and raise the result to the power of n)
where adding and subtracting 1 facilitates taking the limit; similarly, subtracting (which does not change the value of the integral since it integrates, in the principal-value sense, to 0) helps to make subsequent integrals converge. Note that principal value of an integral implies replacing by ; this is tacitly assumed from now on.
Replacing x by ny makes the previous expression into
whose limit is
since and . The last displayed expression is thus the characteristic function of U and (our notation for the corresponding limits of and ).
Note that only the tail behaviour of (8) was relevant in the end.
3. Finding CHF of (1)
Replacing t by (with w positive) results in
where (after the substitution)
Note that the value of is thus always negative (this remains true for the real part of after we make w complex; this becomes consequential later on).
since moving the path of the t integration does not change the integral’s value (there are no singularities between the two paths) and the integrand tends to 0 sufficiently fast within the same strip when w tends to plus or minus infinity.
Therefore, is given by
using the substitution. Note that , since and
This is proved by replacing the path of integration (the real axis) by a clockwise half circle of radius R, centered on the origin and denoted (legitimate, since the integrand has no singularities between the two paths), and then evaluating
The first integral is equal to 0 due to Jordan’s lemma. In the second integral, we have traded the clockwise half circle for counter-clockwise (by changing the integrand’s sign).
Evaluating (15) at two distinct values of w (let us denote them w and w0), dividing the difference by s, and integrating over s from 0 to infinity yields
where is the joint PDF of U and V; note that its support is the whole upper half of the u-v plane.
Replacing s by (note that v cancels out from ) changes the left hand side of (23) to
where is the characteristic function of , the variable whose distribution we seek; the last expression is then equal to the right-hand side of (23).
4. Converting to PDF
Differentiating the resulting equation, i.e. (25) = (23), with respect to w then yields (after a sign reversal)
(to the last integrand, is just a negative constant).
Finally, multiplying both sides of the previous equation by
and integrating over w from 0 to (a notation which implies following the positive imaginary axis) results in
where denotes the corresponding principal value, and (16) is now extended to complex arguments.
Proof. Following  , we write
A further substitution makes it into
Adding (29) and (31), which are identical in terms of their value, and dividing by 2 yields
Finally, introducing results in
assuming that, to evaluate the last integral, we first replace t by (to make it converge), and then take the limit of the answer (Cauchy-type integration).
It is well known that the real part of the left hand side (and, therefore, of the right-hand side) of (28) yields the desired (clearly symmetric) PDF of , say , due to being real (and thus automatically symmetric as well). This then yields
based on (16) and (18). The problem is that there is no simple analytic answer to the last integral; worse yet, its integrand is highly oscillatory, preventing us from direct numerical integration.
In an attempt to break this impasse, we introduce the following substitution
where now follows the complex ray at −45˚.
Unfortunately, even after this substitution, special measures are still needed to facilitate accurate numerical integration of the integral. For one, it becomes necessary to separately deal with the and regions.
5. Case of y2 < 1
When , it is legitimate to rotate the ray of the last integration to the positive real axis, since there are no singularities of the integrand between the old and the new path, and the function decreases sufficiently fast towards infinity in that segment. We then get
where is real, and the integrand is “well behaved” (no oscillations). Results of numerical evaluation are displayed in Figure 1 (we are showing only the portion of the graph; visualize a mirror image of this curve when ).
To investigate the nature of this singularity we notice that, for large , the integrand of (37) tends to
Figure 1. Resulting PDF for . Note the logarithmic-type singularity at .
since . Integrating the last function over from 0 to infinity yields
where the second factor is a special case of the incomplete gamma function defined by
and is Euler’s gamma. This indicates that, by adding to makes the resulting function finite (and thus amendable to a simple, Padé-type approximation, constructed later on).
6. Case of y2 > 1
The situation becomes somehow more difficult when . Now, we can rotate the −45˚ ray of the (36) integration to −90˚ (the negative imaginary axis), since the integrand decreases sufficiently fast towards infinity in this range of directions. The new path contributes 0 to the real part of (36), since the integrand (along the −90˚ ray) is real and is purely imaginary. But moving from −45˚ to −90˚ we have crossed infinitely many singularities located between the two rays; these can be found numerically (a rather tedious procedure) by solving
The location of the first 50 of these (they are all simple poles, slowly converging to the line) are shown in Figure 2.
We can then evaluate (36) using Cauchy integral theorem; note that, at each of these poles (say at ), the derivative of equals to , since
Figure 2. First 50 roots of (41).
This means that each pole (with the exception of the first one—see the next paragraph), contributes
The pole on the imaginary axis needs to be avoided by an infinitesimal half circle, which means that it contributes only one half of the above amount, namely
The last function yields the asymptotic behaviour of as , and provides an excellent approximation (its maximum absolute error is about ) to this PDF when . Unfortunately, to build a numerical solution for the full range of values by adding contributions of sufficiently many of these poles (as done by  ) is rather tedious, as finding hundreds of these poles (needed to reach a good accuracy, especially when approaches 1) is a non-trivial task, and the resulting convergence is quite slow.
Nevertheless we would like to mention that, by adding the contribution of the second pole, namely
where , to (46) results in an equally accurate approximation (its error is less than ) in the tails of .
7. Alternate Solution
We are now going to backtrack and deal directly with (36). Even though its integrand is still highly oscillatory (with a frequency which increases as ), we can mitigate the problem by dividing the range of integration into two segments: first we integrate from 0 to (which does not pose any numerical difficulty), thus getting what we will call the first component of (36), then from to infinity (along the −45˚ ray), getting the second component.
To carry out the latter integration, we first note that, for large , can be expanded in the following manner
Substituting this into the integrand of (36) and further expanding in powers of (up to and including the 6th power) while keeping the two exponential terms fixed (i.e. not expanding in ) yields
which can be -integrated, from to , analytically; then we numerically integrate (over the same line segment) the difference between the integrands of (36) and (49)—since the resulting integrand now approaches zero (as increases) very quickly, the integration range becomes effectively finite ( to already achieves very high accuracy), thus eliminating any troublesome oscillation. Adding the real part of the last two answers (and multiplying by ) then provides the second component of .
This technique yields the graph of Figure 3 (drawn for ; visualize a mirror image of this curve when ):
8. Monte-Carlo Verification
One can always get a good idea about a distribution of any sample statistic (however complicated and inaccessible to analytic treatment) by actually generating a random sample with the required properties by a computer, and using it to compute the desired statistic’s value; this is then repeated as many times as possible, displaying the results in a histogram.
We have done this with our (1), using and generating one million of its random values. Plotting the resulting histogram, together with the theoretical asymptotic PDF of the last three sections, yields Figure 4; visual comparison clearly indicates a good agreement between the two answers. The tiny discrepancy still discernible (mainly in the range) is due to the fact that is not large enough to have reached the limit yet (trying to make the sample size substantially higher would run up against our computer’s capacity).
Note that, to be more economical, we have folded the negative and positive parts of the distribution into a single graph, effectively plotting the PDF of the absolute value of (1).
9. Accurate Approximation
First we have to remember the we already have an excellent approximation to in the region, given by the sum of (46) and (47); now we have to find a similarly accurate way of dealing with .
Figure 3. Resulting PDF when .
Figure 4. Comparing theoretical and empirical PDF’s.
We already know that, by adding
to removes the corresponding singularity when ; we can easily verify that the same is true when , since one of the terms in the analytic result of the (49) integration (once we take its real part and multiply by ) equals to
whose singular part is given by (50). This proves that
is singularity-free for all values of y; unfortunately, that still does not make it sufficiently “smooth”, as the next paragraph indicates.
There is yet another subtle issue causing great difficulty when trying to build an approximate formula for (53): the function has a small (hardly noticeable) kink (a discontinuous second derivative) at ; this is not an artifact of the new technique—the same (rather surprising) phenomenon can be confirmed by the old technique of adding residues. Luckily, we can identify yet another term of the (49) contribution responsible for this discontinuity, and remove it by further subtracting the offending term from (53), thus getting
This function is (finally!) sufficiently (even though not perfectly—tiny issues remain with higher derivatives at etc.) smooth to fit it (in the range) by a ratio of two polynomials (by minimizing the total error squared). This results in the following approximation
Its absolute error never exceeds , which is more than sufficient for any practical application. The form of this expression and the individual powers of y have been chosen somehow arbitrarily; we do not claim that our choices are optimal.
Do not forget that (55) is approximating (54); a corresponding adjustment has to be made to convert it into an approximation for .
The aim of this article was to demonstrate that finding the distribution of a relatively simple sample statistic requires a skillful use of characteristic functions and a whole gamut of sophisticated mathematical techniques, including real and complex analysis, Fourier transform, and curve fitting. We hope that students of Statistics can benefit from the ingenuity of the authors of the original derivation of this distribution as presented in  , and from some extra details included in this article; the latter includes a different numerical approach to building the resulting PDF, expressing it in the form of an accurate Padė-type approximation (discovering an interesting discontinuity in the process), and verifying the answer by Monte Carlo simulation.