A Note on “Limit Distributions of Self-Normalized Sums” Using Cauchy-Generated Samples

In this case study, we would like to illustrate the utility of characteristic functions, using an example of a sample statistic defined for samples from Cauchy distribution. The derivation of the corresponding asymptotic probability density function is based on [1], elaborating and expanding the individual steps of their presentation, and including a small extension; our reason for such a plagiarism is to make the technique, its mathematical tools and ingenious ar-guments available to the widest possible audience.


Introduction
The problem of finding the distribution of self-normalized sum of a Cauchygenerated sample has been considered since at least 1969 (see, for example, [2]), but solved (by proposing a concrete numerical algorithm) only in 1973 by [1], 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, [3].
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 n → ∞ limit of the probability density function (PDF) of How to cite this paper: Vrbik, J. (2019) A Note on "Limit Distributions of Self-Normalized Sums" Using Cauchy-Generated where 1 2 , , , n X X X 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 [1]; 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 ( ) f x stands for the X's PDF,  for the corresponding expected value, and i is the purely imaginary unit; note that the real part of ( ) x s ϕ is always an even (an alternate name is symmetric) function of s, while the purely imaginary part is odd. Also note that when ( ) f x is even, the corresponding CHF must be real (its purely imaginary part becomes zero because the integrand is odd

( )
, f x y 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 X Y + is given by , further implying that the CHF of a sample mean of n independent i X is given by When two random variables are combined into one (e.g. defining X U Y = ), 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 [4].

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 2 X (seen as distinct random variables) is then given This implies that the joint CHF of def 1 and of is (replace s by s n and t by 2 t n in (9), and raise the result to the power of n) where adding and subtracting 1 facilitates taking the n → ∞ limit; similarly, subtracting ixs n (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 dx The last displayed expression is thus the characteristic function of U and 2 V (our notation for the corresponding limits of n U and 2 n V ).
Note that only the tail behaviour of (8) was relevant in the end.

Finding CHF of (1)
Replacing t by 2 is w (with w positive) results in where (after the ys z = substitution) Note that the value of ( ) w ψ is thus always negative (this remains true for the real part of ( ) w ψ after we make w complex; this becomes consequential later on). Proof.
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, ( ) w ψ is given by 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 eva- 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 w 0 ), dividing the difference by s, and integrating over s from 0 to infinity yields f u v 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
, the variable whose distribution we seek; the last expression is then equal to the right-hand side of (23).

Converting to PDF
Differentiating the resulting equation, i.e. (25) = (23), with respect to w then yields (after a sign reversal) Finally, multiplying both sides of the previous equation by and integrating over w from 0 to i ⋅ ∞ (a notation which implies following the positive imaginary axis) results in after introducing Adding (29) and (31), which are identical in terms of their value, and dividing by 2 yields assuming that, to evaluate the last integral, we first replace t by 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 U V , say ( ) f y , due to ( ) 2 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 2 1 y < and 2 1 y > regions.

Case of y 2 < 1
When 2 1 y < , 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  To investigate the nature of this singularity we notice that, for large τ , the integrand of (37) tends to where the second factor is a special case of the incomplete gamma function de-

Case of y 2 > 1
The situation becomes somehow more difficult when 2 1 y > . 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 dτ 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 Im Re τ τ = − line) are shown in Figure 2. We can then evaluate (36) using Cauchy integral theorem; note that, at each of these poles (say at p τ ), the τ derivative of ( )   ( ) This means that each pole (with the exception of the first one-see the next paragraph), contributes to (36).
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, The last function yields the asymptotic behaviour of ( ) f y as 2 y → ∞ , and provides an excellent approximation (its maximum absolute error is about 6 2 10 − × ) to this PDF when 1.8 y > . Unfortunately, to build a numerical solution for the full range of 2 1 y > values by adding contributions of sufficiently many of these poles (as done by [1]) is rather tedious, as finding hundreds of these poles (needed to reach a good accuracy, especially when 2 y 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 Applied Mathematics Im exp 2 where 2 1.84906 3.45208i τ = − − , to (46) results in an equally accurate approximation (its error is less than 6 2 10 − × ) in the 1.7 y > tails of ( ) f y .

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 To carry out the latter integration, we first note that, for large τ , ( ) Q τ can be expanded in the following manner ( ) Substituting this into the integrand of (36) and further expanding in powers of 1 τ which can be τ -integrated, from 1 i − to ( ) 1 i ∞ − , 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 (1 i

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 clearly indicates a good agreement between the two answers. The tiny discrepancy still discernible (mainly in the 0 1 y < < range) is due to the fact that 100 n = is not large enough to have reached the n → ∞ 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).

Accurate Approximation
First we have to remember the we already have an excellent approximation to ( ) f y in the 1.7 y > region, given by the sum of (46) and (47); now we have to find a similarly accurate way of dealing with 1.7 y ≤ .   to ( ) f y removes the corresponding singularity when 2 1 y ≤ ; we can easily verify that the same is true when 2 1 y > , since one of the terms in the analytic result of the (49) integration (once we take its real part and multiply by 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 2  Its absolute error never exceeds 6 6 10 − × , 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 ( ) f y .

Conclusion
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 [1], 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.

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.