Efficient computation of the cumulative distribution function of a linear mixture of independent random variables

TL;DR

This paper proves the correctness of an improved algorithm for efficiently computing the cumulative distribution function of a linear mixture of independent random variables, with applications in bootstrap methods and hypothesis testing.

Contribution

It provides a proof of convergence for a variant of an existing algorithm, ensuring accurate approximation of the distribution function from finite samples.

Findings

The algorithm converges to the true distribution as resolution increases.

It efficiently computes the complete distribution function from sample data.

Applications include bootstrap estimates and hypothesis testing in linear models.

Abstract

For a variant of the algorithm in [Pit19] (arXiv:1903.10816) to compute the approximate density or distribution function of a linear mixture of independent random variables known by a finite sample, it is presented a proof of the functional correctness, i.e. the convergence of the computed distribution function towards the true distribution function (given the observations) as the algorithm resolution is increased to infinity. The algorithm (like its predecessor version) bears elements which are closely related to early known methods for numerical inversion of the characteristic function of a probability distribution, however here efficiently computes the complete distribution function. Possible applications are in computing the distribution of the bootstrap estimate in any linear bootstrap method (e.g. in the block bootstrap for the mean as parameter of interest, or residual bootstrap in linear regression with fixed design), or in elementary analysis-of-variance hypothesis testing.

Equations

$$Z:=\sum _{j=1}^m{a}_j{X}^{[j]}$$

$$G_{a_j{X}^{[j]}}(\nu )$$

$$G(\nu )$$

$$\tilde{h}(x)$$

$$\tilde{f}(x)$$

$${\hat{f}}_Z(x)$$

$$F_Z=\sum _{z\in \{Z\}}{p}_z\cdot 1_{z\le x}.$$

$$\bar{F}(x):=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{e^{2\pi i\cdot x\nu }-e^{2\pi i\cdot x_0\nu }}{i\nu }G(\nu )\text{d}\nu ,$$

$$H(x)$$

$$F_Z(x)=\bar{F}(x)$$

$${\tilde{F}}_{i^{\prime }}$$

$$=\frac{1}{N}\sum _{k=-N+1}^{N-1}G(k\Delta \nu )\sum _{i^{\prime \prime }=0}^{i^{\prime }-i_{min}-1}\exp (2\pi i\cdot \frac{(i_{min}+i^{\prime \prime })k}{N})$$

$$\tilde{f}=\frac{T}{N}\cdot \sum _{z\in \{Z\}}{p}_z\cdot R_{N,T}((\cdot )-z)$$

$$R_{N,T}(x):=\frac{1}{T}\cdot \frac{\sin (2\pi \frac{2N-1}{2T}\cdot x)}{\sin (2\pi \frac{1}{2T}\cdot x)}$$

$$\tilde{f}(x)=\frac{T}{N}\cdot \left(f_Z\ast R_{N,T}\right)(x)$$

$$f_Z=\sum _{z\in \{Z\}}{p}_z\cdot {\delta }_z$$

$$\tilde{f}(x)$$

$$=\frac{T}{N}\sum _{z\in \{Z\}}{p}_z\cdot \frac{1}{T}\sum _{k=-N+1}^{N-1}\exp (2\pi i\cdot (x-z)k/T)$$

$$1\stackrel{!}{=}{\int }_{\kappa \cdot z_{min}}^{\kappa \cdot z_{min}+T}\hat{f}(x)\text{d}x={\int }_{\kappa \cdot z_{min}}^{\kappa \cdot z_{min}+T}c\cdot N\cdot \tilde{f}(x)\text{d}x\approx c\cdot N\cdot \sum _{i=0}^{N-1}{\tilde{f}}_i\cdot \frac{T}{N},$$

$$\hat{f}(x):=\frac{N}{T}\tilde{f}(x)\cdot 1_I(x)$$

$$\tilde{F}(x):={\int }_{-\infty }^x\frac{N}{T}\cdot \tilde{f}(\xi )\cdot 1_I(\xi )\text{d}\xi .$$

$$\tilde{F}(x)\to F_Z(x)\forall x\in \mathbb{R}\backslash A_d$$

$$\underset{N\to \infty }{\lim }\frac{N}{T}\cdot \tilde{f}$$

$$\frac{N}{T}\cdot \tilde{f}\cdot 1_I\to (f_Z\ast \text{LLI})\cdot 1_I=f_Z.$$

$$\underset{N\to \infty }{\lim }$$

$$={\int }_{I_x}(f_Z\ast \text{LLI})=F_Z(x)$$

$${\tilde{F}}_i-\tilde{F}(iT/N)\to 0\forall i=0\dots N-1$$

$$\tilde{F}(x)$$

$$=\frac{x-x_0}{T}+\frac{1}{T}\cdot \sum _{0<|k|\le N-1}{g}_k{\int }_{x_0}^x\exp (2\pi i\cdot \xi k\Delta \nu )\text{d}\xi$$

$$=\frac{x-x_0}{T}+\sum _{0<|k|\le N-1}{g}_k\frac{\exp (2\pi i\cdot xk\Delta \nu )-\exp (2\pi i\cdot x_{0}k\Delta \nu )}{2\pi i\cdot k}.$$

$$\tilde{F}(iT/N)$$

$$=\frac{i-i_{min}}{N}+\frac{1}{N}\cdot \sum _{0<|k|\le N-1}{g}_k\cdot \frac{e^{2\pi i\cdot ik/N}-e^{2\pi i\cdot i_{min}k/N}}{2\pi i\cdot k/N}$$

$$=\frac{i-i_{min}}{N}+2\text{Re}\left(\frac{1}{N}\cdot \sum _{k=1}^{N-1}{g}_k\cdot \frac{e^{2\pi i\cdot ik/N}-e^{2\pi i\cdot i_{min}k/N}}{2\pi i\cdot k/N}\right).$$

$${\tilde{F}}_i$$

$$=\frac{i-i_{min}}{N}+\frac{1}{N}\sum _{0<|k|\le N-1}{g}_k\cdot \frac{e^{2\pi i\cdot \frac{ik}{N}}-e^{2\pi i\cdot \frac{i_{min}k}{N}}}{\exp (2\pi i\cdot k/N)-1}$$

$$=\frac{i-i_{min}}{N}+2\text{Re}\left(\frac{1}{N}\sum _{k=1}^{N-1}{g}_k\cdot \frac{e^{2\pi i\cdot \frac{ik}{N}}-e^{2\pi i\cdot \frac{i_{min}k}{N}}}{\exp (2\pi i\cdot k/N)-1}\right).$$

$${\widetilde{\stackrel{~}{F}}}_i=\tilde{F}(iT/N).$$

Full text

Keywords: characteristic function inversion, linear independent mixture, variance reduction

AMS 2010 Subject Classification: 62G07, 62G09

1 Introduction

Bootstrapping is a resampling technique employed to estimate a parameter $\theta$ of a distribution (most often variances, in order to systematically construct confidence intervals) in the face of no a-priori knowledge about the true distribution from which the available data points $X_{i}$, $i=1\dots n$, (in this text real-valued) are deemed to be originating from. In particular, bootstrapping is capable of extending hypothesis testing beyond the requirement that a parametric form of the distribution of the involved noise is known.

The conventional way of using bootstrapping in practice is to evaluate a certain estimator $T_{n}^{*}$ on multiple ”synthetic” data sets (=replicates) which are generated (in the basic type) by randomly sampling data points from the orginal data set, with the aim of letting them mimick this original data set with regards to its distributional characteristics. The so derived values will, under suitable assumptions placed on the bootstrap method (i.e. on estimator $T_{n}^{*}$ and the resampling distribution) as well as on the parameter of interest, mimick the distribution of an estimator $T_{n}$ which is a centered and scaled derivate of the estimator $\hat{\theta}$ of the actually targetted parameter $\theta$. (=”consistency”) The method thus allows, under suitable assumptions, to find the distribution of $\hat{\theta}$ (when applied to the whole population) and thus the ultimately targetted confidence regions.

In [Pit19], an algorithm was outlined which allows to circumvent the above described ”random element” in conventional bootstrap method usage and arrive at an approximation of the distribution of $T_{n}^{*}$ (conditional on the available data points $X_{i}$) by directly (deterministically) concluding from the $X_{i}$, provided that the bootstrap method used is linear in the sense defined in [Pit19]. In the present text, the aim is to prove rigorous statements on the error of the approximation, and the focus will be on an in $N$ asymptotic result. The main statement on this is contained in section 4.2, as well as in the appendix.

Besides facilitating practical realization of the bootstrapping procedure for linear bootstrap methods, the computation of the distribution of a linear mixture of random variables potentially has applications related to other estimation problems. In the wider context of parameter estimation in econometric models at small sample sizes, [Phi82] surveys methods for estimating the distribution of estimators in simultaneous equations models. [Dav80] has error bounds on the (for this field relevant) computation of the distribution of quadratic forms of multiple independent normally distributed random variables, which may be written as linear combinations of independent $\chi^{2}$ variables, i.e. which are in a form amenable to the here examined algorithm. In [ACW99], methods are developed for inverting transforms of probability distributions, and applied to obtain tail probabilities in queuing systems. (Compared to their work, which is geared towards obtaining single probability values with high accuracy and typically uses analytical expressions of the characteristic functions of the underlying pdfs, the algorithm here computes the whole distribution function (sampled equidistantly).)

The field of most direct application could initially appear to be the hypothesis testing for the comparison of two populations’ means, which occurs throughout social and life sciences, see for example [MH03], [HWC*+*16]. However, the recognized established non-parametric tests (Mann-Whitney test, resp. rank-based Kruskal-Wallis test in the more-than-two-groups case) will be preferable due to the more rigorous statements derivable at equal assumptions. Such test scenarios may be regarded as analysis-of-variance with discretely levelled factors. For confidence intervals in factor models with continuously valued factors, applicability of the here examined algorithm remains an open question.

The text proceeds as follows: In the next section, the problem to be tackled is restated, and elementary definitions and relations given. In section 3, elementary properties relevant to the inversion of the characteristic functions are stated. In section 4, an alternative explicit expression of the algorithm output is derived and the main convergence result proven with it.

2 Setting and common definitions

Let $X_{i}$, $i=1\dots n$, be real-valued observations deemed to be realizations from some (not further considered) random variable, and let $\hat{F}_{n}$ be the associated empirical distribution function, i.e. $\hat{F}_{n}(x)=n^{-1}\cdot\sum_{i=1}^{n}1_{X_{i}\leq x}$. Let the random variable $X$ be defined as distributed according to $\hat{F}_{n}$. Let $X^{[j]}$, $j=1\dots m$, be independent random variables distributed as $X$. Let $a_{j}\in\mathbb{R}$, $j=1\dots m$. Define

[TABLE]

and denote by $F_{Z}$ its cumulative distribution function. In [Pit19], it was stated that the algorithm presented there, here called Algorithm 1, computes an approximation to the density $f_{Z}$ of $Z$. With $G$ and $g_{k}$ defined as222It is ”$i$” after ”$2\pi$” or in denominator the imaginary unit, otherwise usually acts as index variable.

[TABLE]

i.e. $g_{k}$ as computed in Algorithm 1, set for $x\in\mathbb{R}$ and $N\in\mathbb{N}$

[TABLE]

It is $g_{0}=1$ and $\left\lvert g_{k}\right\rvert\leq 1$; because of $g_{0}=1$ and $g_{-k}=g_{k}^{*}$, one easily finds $\tilde{f}(x)=2\cdot\text{Re}(\tilde{h}(x))-1/N$. For an $i\in\{0,\dots,N-1\}$, it is $\tilde{f}_{i}:=\tilde{f}(i\tfrac{T}{N})$ (called $f_{i}$ in equation (5) in [Pit19]). Set $\tilde{h}_{i}:=\tilde{h}(i\tfrac{T}{N})$. One comfortably proves $\sum_{i=0}^{N-1}\tilde{h}_{i}=1$ and therefore also $\sum_{i=0}^{N-1}\tilde{f}_{i}=1$.

An alternative representation of $F_{Z}$ can be stated by recognizing that $Z$ takes finitely many values. Denote by $\{Z\}$ the set of those values, and let $p_{z}$ the probability mass of value $z\in\{Z\}$ (as implied by $\hat{F}_{n}$). Then

[TABLE]

3 Towards convergence proof: Elementary properties of $G$

In this section, the proof of the convergence of the algorithm result to the desired distribution $F_{Z}$ is prepared by recollecting some elementary properties for the inversion of a characteristic function and stating relevant definitions. ”Convergence” and ”asymptotic” here refers to the behaviour as $N\rightarrow\infty$, while (on the contrary) the $n$, i.e. the number of data points, remains fixed. Insofar, it is examined here the necessary computational resources to be expended to achieve a sufficiently accurate result, ideally independent of the given fixed input sample size $n$.

In prospect of wanting to apply the knowledge of the distribution to deriving confidence intervals and rejection probabilities, it is the error in estimating $F_{Z}(z)$ for each $z\in\mathbb{R}$ that is of interest.

Let $G(\nu)$ be as in eqn. (3). Remark 3 below asserts that $G$ equals –up to argument-side scaling– the characteristic function of the distribution of $Z$. Since the characteristic function determines the distribution uniquely, it is $G(\nu)$ representing this distribution exactly. Applying then an inversion formula, the earliest version of which appears to have been stated by Lévy and a modified version of which is stated below as Lemma 3, therefore yields the true cumulative distribution function $F_{Z}$ of $Z$.

Remark 1:

Let $X_{i},i=1\dots n,$ be real-valued observations, and $\hat{F}_{n}$ the associated empirical distribution function. Let random variable $Z$ be as in equation (1), again with the contained random variables $X^{[j]}$ independent and distributed according to $\hat{F}_{n}$. Let $G$ be defined as in (3). Then the characteristic function of $Z$ (conditional on $X_{1},\dots,X_{n}$), defined as $t\mapsto\text{E}(\exp(itZ))$, is given by $t\mapsto G(-t/(2\pi))$.

The statement of the remark is based on the well-known argument that by independence of the $X^{[j]}$, the expectation separates into a product of expectations. Since the $X^{[j]}$ are discretely distributed, each of the expectations can be written as in (2). ∎

The following lemma is stated here for completeness, but will not be made use of directly. The expressions derive from the well-known inversion formulae for the characteristic function originally appearing in [Lév25], by simple variable substitution and noting that $F_{Z}(x_{0})=0$.

Lemma 1:

Let $X_{i}$, $\hat{F}_{n}$, $Z$ and $G$ as in Lemma 1. Let $x_{0}<\min(Z)$. Then at every point $x$ of continuity of $F_{Z}$, it is $F_{Z}(x)=\bar{F}(x)$ with

[TABLE]

and because $G(-\nu)=G^{*}(\nu)$, at the same points also $F_{Z}(x)=2\cdot\text{Re}(H(x))$ with

[TABLE]

Here, since $\left\lvert G(\cdot)\right\rvert$ is bounded and the fraction in the integrand of (9) tends to $2\pi(x-x_{0})$ as $\nu\rightarrow 0$, the integrand is continuous and bounded everywhere and the integral in (9) can well be evaluated as Riemann integral around $\nu=0$.

∎

One route of analysis at this point would proceed by using Lemma 3 to establish a link between an estimate of the distribution function (defined in the second next paragraph) and the true $F_{Z}$ (as given by its integral representation above). For the purpose of this text however, a more accessible and possibly more illuminating route is preferred.

The following definitions are used throughout the remainder of the text. By $\tilde{F}_{i}$ (see equation (11) below) it is denoted the ”cyclical” accumulative sum of the $\tilde{f}_{i}$ of Algorithm 1. It maps ”$\!\!\!\!\mod N$” into the integer range $\{0,\dots,N-1\}$.

Set $T_{Z}:=\max Z-\min Z$. Let $T>T_{Z}$, $\kappa:=T/T_{Z}$, assume $z_{min}<0<z_{max}$ and let $i_{min}=\left\lfloor\tfrac{N}{T}\cdot\kappa\cdot z_{min}\right\rfloor$ and then $x_{0}=i_{min}T/N$. For $i^{\prime}=i_{min}\dots(i_{min}+N-1)$ set

[TABLE]

4 Results at finite $N$ and asymptotic result

The route of analysis pursued in this section aims to represent the algorithm output in terms of desired or known quantities. The argument has similarity to ones known in the Fourier theory of sampling.

Theorem 1:

Let $X_{i}$, $\hat{F}_{n}$, $Z$, $G$ and $\tilde{f}$ be given as in section 2. Then

[TABLE]

where

[TABLE]

with the expression on the right deemed continuously continued at the zeros of the denominator.

Remark 2: The statement of the theorem is also written as

[TABLE]

where ”$*$” signifies the convolution of two tempered distributions (i.e. in the sense of a generalized function, see appendix A) , and

[TABLE]

is the tempered distribution representing the probability density belonging to $F_{Z}$. (The $\delta_{z}$ are delta distributions with mass at $z$.)

Proof of the theorem:

Utilizing the representation of $F_{Z}$, it is the Fourier transformation $G$ of $Z$ written as $G(\nu)=\text{E}(\exp(-2\pi i\cdot Z\nu))=\sum_{z\in\{Z\}}p_{z}\cdot\exp(-2\pi i\cdot z\nu)$. Substituting this in the algorithm output (5) yields

[TABLE]

Identifying the finite geometric series in the sum over $k$ and applying an exponential factor at it yields the above stated expression for $R_{N,T}$. (The algorithm output is obtained by evaluating $\tilde{f}$ at the places $i\tfrac{T}{N}$, $i=0\dots N-1$.) ∎

The $R_{N,T}$ has properties which will allow to derive useful characteristics of the algorithm output. It is $\lim_{x\rightarrow 0+kT}R_{N,T}(x)=\tfrac{2N-1}{T}$ for $k\in\mathbb{Z}$. It is $R_{N,T}$ periodic with period $T$, and $\int_{0}^{T}R_{N,T}(x)\,\text{d}x=1$. As $N\rightarrow\infty$, have R_{N,T}\rightarrow\raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T} in the space $\mathcal{S}^{\prime}$ of tempered distributions (see appendix A), where \raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T} denotes the Dirac comb333Informal definition: \raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T}:=\sum_{k=-\infty}^{\infty}\delta((\cdot)-kT). with interval $T$. (The envelope of $R_{N,T}$, as $N$ increases, remains constant, but the increasing number of zero crossings means that in intervals not overlapping with $T\cdot\mathbb{Z}$, the positive and negative contributions cancel.)

4.1 Motivating the proposed expression for the density estimate

According to Theorem 4, the quantity $\tilde{f}$ can be regarded as proxy for a density estimate, and (according to the remark subsequent to the common definitions) itself fulfills a normalizing constraint $\sum_{i=0}^{N-1}\tilde{f}_{i}=1$. Because this constraint implies that the $\tilde{f}_{i}$ scale to zero like $N^{-1}$ as $N\rightarrow\infty$, a to-be-defined density estimate $\hat{f}$ must reasonably be of form $c\cdot N\cdot\tilde{f}$ in order to asymptotically fulfill the normalization. The $c$ is found by noting that (or rather: aiming for)

[TABLE]

where the $\approx$ sign appeals to the numerical integration of the $\tilde{f}$ integral using $N$ equidistant samples at $i\tfrac{T}{N}$. Thus reasonably $c=1/T$.

Thus, if we here and henceforth set $I=[\kappa\cdot z_{min},\kappa\cdot z_{max}]$, then the reasonable estimate for a smooth approximation to the density of $Z$ is

[TABLE]

4.2 Asymptotics

It shall now be proven, using the representation as periodic superposition of the smoothed true density, that the integrated Algorithm 1 output converges to the true cumulative distribution function as $N\rightarrow\infty$. With ”integrated” initially is meant the appropriate summation of the $\tilde{f}_{i}$ (regardable as scaled density estimate according to the previous proposition).

Theorem 2:

Let $X_{i},a_{j},Z,f_{Z},F_{Z},T_{Z}$ be given as in Theorem 4, let $T$ be chosen $T>T_{Z}$, and $N$ be chosen. Let $\tilde{f}$ be defined as in equation (5) (equalling the expression (15) below Theorem 4). Call $I:=[\kappa\cdot z_{min},\kappa\cdot z_{max}]$ and $I_{x}:=I\cap[\kappa\cdot z_{min},x]$ (with $I_{x}=\emptyset$ for $x<\kappa\cdot z_{min}$), and set

[TABLE]

Let $A_{d}$ be the set of points of discontinuity of $F_{Z}$. Then it holds as $N\rightarrow\infty$:

[TABLE]

Remark: The theorem (as the previous ones) considers the $X_{i}$ as given. Consequently also $F_{Z}$ is given as deterministic quantity. The convergence therefore is rightfully meant as a deterministic one.

Proof of (i): As before, $\mathcal{S}$ shall denote the (Schwartz) space of rapidly decreasing functions, and $\mathcal{S}^{\prime}$ the associated space of tempered distributions. It was already stated that R_{N,T}\rightarrow\raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T} in the distributional limit sense as $N\rightarrow\infty$. Further, the convolution operation (among distributions) is continuous in the sense here needed444See results from appendix A. It follows

[TABLE]

Since for $T>T_{Z}$ the expression f_{Z}*\raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T} denotes the periodic repetition of $f_{Z}$, it follows

[TABLE]

Next, note that $F_{Z}(x)=\int_{-\infty}^{x}f_{Z}$ for all $x\not\in A_{d}$, where the integral of the distribution $f_{Z}$ over $[a,b]$ is defined via evaluation at a mollified $1_{[a,b]}$ (see appendix A). For any sequence of $\mathcal{S}$ functions $g_{n^{\prime}}$ which converge (embedded in $\mathcal{S}^{\prime}$) towards a distribution $G$, it is then $\int_{a}^{b}(g_{n^{\prime}})\rightarrow\int_{a}^{b}G$ as ${n^{\prime}}\rightarrow\infty$. Consequently here also

[TABLE]

for all $x\in\mathbb{R}\backslash A_{d}$. This completes the proof of i). ∎

It is thus seen that the integral behind $\tilde{F}(\cdot)$ tends to the desired function $F_{Z}$ at all places of continuity of $F_{Z}$. One would now want to continue by proving a statement along

[TABLE]

under the same assumptions as used for i), and where $\tilde{F}_{i}$ is the sum of the $\tilde{f}_{i}$ as in equation (11). The viewpoint behind this claim is that, in schematic words, $\tilde{F}_{i}$ is the sum of the $\tilde{f}_{i}=\tilde{f}(i\tfrac{T}{N})$ over suitable interval, while $\tilde{F}(\cdot)$ is the integral of $\tilde{f}(\cdot)$ over the same interval. Thus $\tilde{F}_{i}$ roughly is the numerical evaluation (by equidistant sampling and rectangular rule) of the integral behind $\tilde{F}(\cdot)$. A problem here occurs since as $N$ increases and the integration partitioning becomes finer, also the integrand changes. It is then possible that the evaluations $\tilde{f}_{i}$ of $\tilde{f}$ all are errorenous, with the errors accumulating in the sum representing $\tilde{F}_{i}$. This will occur particularly if the possible values of $Z$, and thus the placement of the $R_{N,T}(\cdot)$, exhibit a regularity. (Unrelated to the regularity in $f_{Z}$, the sampling of a smoothed $f_{Z}$ can be expected to produce errorenous results for certain choices of $N$ in a way that exhibits oscillation as $N$ varies.)

To naturally avoid this, it is more suitable to sample $\tilde{F}(\cdot)$. With $x_{0}:=\kappa\cdot z_{min}$, have

[TABLE]

Evaluated at $x=iT/N$, and using $x_{0}=i_{min}T/N$, obtain

[TABLE]

This looks remarkably similar to the expression already available by the summation of the $\tilde{f}_{i}$. It is

[TABLE]

One concludes that the desired values of $\tilde{F}(iT/N)$ can be obtained by computing the $\tilde{f}_{i}$ as in the previous way, but with the $g_{k}$ suitably modified by a factor.

Concretely, denoting by $\tilde{\tilde{F}}_{i}$ the expression constructed from $\tilde{F}_{i}$ by replacing $G(k\Delta\nu)$ in (11) with $G(k\Delta)\cdot(\exp(2\pi i\cdot k/N)-1)/(2\pi i\cdot k/N)$, using the limit ”$k\rightarrow 0$” of the fraction at $k=0$, one obtains

[TABLE]

This is recorded in the lemma below.

The modified version of the Algorithm 1 is henceforth referred to as ”Algorithm 2”.

Lemma 2:

Let $X_{i},a_{j},Z,G,f_{Z},F_{Z},T_{Z}$ be given as in Theorem 4.2, let $T$ be chosen as $T>T_{Z}$. Let $\tilde{f}$ be defined as in (5), and $\tilde{F}$ as in (21). Let $\tilde{\tilde{F}}_{i}$ be as defined before equation (32). Then

[TABLE]

Proof: As argued before the statement of the lemma.∎

Using i) of Theorem 4.2 and the above ii’), one concludes that for all $i=0\dots N-1$ with $x=i\tfrac{T}{N}\not\in A_{d}$, it holds $\tilde{\tilde{F}}_{i}\rightarrow F_{Z}(i\tfrac{T}{N})$ as $N\rightarrow\infty$. Clearly, using $N$ sufficiently large, it will suffice to evaluate $\tilde{F}$ at those discrete places for which the theorem does provide the convergence assertion. In practice and when not targetting specific assertions on the error, a value of $N=1000$, in some applications $N=10000$, appears sufficient to compute for example quantiles of the distribution of $Z$. A derivation of a rigorous bound on the error at finite $N$ is found in the appendix B.

5 Conclusion

The convergence of a variant of the algorithm in [Pit19] (see modification derived in section 4.2) has been proven and a bound on the absolute value of error stated in dependence on the algorithm resolution $N$. Future research could target further improvement of the convergence behaviour, in particular by more effectively using the knowledge about the periodicity in the error component. Moreover, details of the application of the algorithm to the areas mentioned in the introduction could be examined.

Appendix A Functional-analytic background for proof of Thm. 4.2

The objects in equation (15) and in the proof of Theorem 2 part (i) are deemed elements of $\mathcal{S}^{\prime}$, the set of tempered distributions [RS72a], i.e. continuous linear functionals $\mathcal{S}\rightarrow\mathbb{C}$. Here $\mathcal{S}$ is the space of rapidly decreasing $C^{\infty}$ functions $\mathbb{R}\rightarrow\mathbb{C}$, endowed with the family of semi-norms $\lVert\phi\rVert_{\alpha,\beta}:=\max_{\mathbb{R}}\left\lvert x^{\beta}D^{\alpha}\phi\right\rvert$ (see [RS72a]). Even though the in Theorem 4 mentioned $R_{N,T}(\cdot)$ is not in $\mathcal{S}$, regarding $R_{N,T}(\cdot)$ as tempered distribution allows

[TABLE]

with converging integral on the right-hand side, so $\mathcal{F}[R_{N,T}]$ is a well-defined functional on $\mathcal{S}$, and is (because of continuity of $\mathcal{F}$ in $\mathcal{S}$, see [RS72b]) continuous. Thus $\mathcal{F}[R_{N,T}]\in\mathcal{S}^{\prime}$, and moreover its limit in $\mathcal{S}^{\prime}$ is well-defined. (In the following text, reference to ”tempered” will be dropped, even though meant.)

For a distribution $F$ and a $g\in\mathcal{S}$, it is defined $(F*g)(\phi):=(F)((g\circ m)*\phi)$, see [RS72b]. For $F\in\mathcal{S}^{\prime}$, it will be $F*\delta_{a}$ be deemed defined via approximating function sequence $\bar{\delta}_{a,1/s}$, knowing $\bar{\delta}_{a,1/s}\rightarrow\delta_{a}$ as $s\rightarrow\infty$. Then easily confirmed: $F*\delta_{0}=F$, and $f(\cdot)*\delta_{a}=f((\cdot)+a)$, and desirable properties (e.g. the convolution/product theorem) carry over via continuity.

More generally, for any distribution $G$ to which a sequence of $\mathcal{S}$ functions $g_{n}$ converges (in $\mathcal{S}^{\prime}$), can define $(F*G):=\lim_{n\rightarrow\infty}(F*(g_{n}))$. Then $G\mapsto(F*G)$ is continuous for those sequences, i.e. from $\lim_{n\rightarrow\infty}(g_{n})=G$ can conclude $\lim_{n\rightarrow\infty}(F*(g_{n}))=(F*\lim_{n\rightarrow\infty}(g_{n}))=(F*G)$.

The integral of a distribution over a bounded interval is defined as follows: for $a<b$, let $\bar{J}_{[a,b],\epsilon}$ be the function $1_{[a,b]}(\cdot)*\bar{\delta}_{0,\epsilon}$ with $\bar{\delta}_{0,\epsilon}$ as in proof of Theorem 4. It is $\bar{J}_{[a,b],\epsilon}\in\mathcal{S}$. Then set for $F\in\mathcal{S}^{\prime}$, if the limit exists and is independent of the shape of the mollifier choice $\bar{\delta}_{0,\epsilon}$,

[TABLE]

Clearly, for $F=(f)$ i.e. a distribution associated to an $L^{1}_{loc}$ function $f$, it is $\int_{a}^{b}F$ existing and equal to the ordinary Lebesgue integral of $f$ over $[a,b]$.

A.1 Convergence of $R_{N,T}$

In section 4, it was stated without proof that R_{N,T}\rightarrow\raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T} as $N\rightarrow\infty$. This is seen as follows: It was $R_{N,T}(x)=T^{-1}\cdot\sum_{k=-N+1}^{N-1}\exp(2\pi i\cdot k\Delta\nu x)$, therefore get (as distributional equation)

[TABLE]

where the last equation uses a well-known result for the Dirac comb \raisebox{-0.56905pt}{\text{\mathsf{L!L!I}}}_{T}. 555The definition of the Fourier transformation $\mathcal{F}$ on distributions underlying this equation is rooted in the definition of the Fourier transformation on the Schwartz space $\mathcal{S}$, using, for $v\in\mathcal{S}$,

$\displaystyle\mathcal{F}[v](\nu)$ $\displaystyle:=\int_{\mathbb{R}}v(x)e^{-2\pi i\cdot x\nu}\,\text{d}x,\text{ and its inverse}$

(39)

$\displaystyle\mathcal{F}^{-1}[\hat{v}](x)$ $\displaystyle:=\int_{\mathbb{R}}\hat{v}(\nu)e^{2\pi i\cdot x\nu}\,\text{d}\nu.$

(40)

Appendix B Rate of uniform convergence of $\tilde{F}(\cdot)$

B.1 Overview

The section gathers the elementary results leading up to the result

[TABLE]

It shall denote in this appendix $\{Z\}$ the multiset of all possible values of $Z$ (each attained with probability $N_{Z}^{-1}$), and $N_{Z}$ its cardinality. It is then $f_{Z}$ simply $f_{Z}=N_{Z}^{-1}\cdot\sum_{z\in\{Z\}}\delta_{z}$. Let $T>T_{Z}$, $\kappa=T/T_{Z}$, $I=[\kappa\cdot z_{min},\kappa\cdot z_{max}]$ as before, and $x_{0}=\kappa\cdot z_{min}$. By equation (21), it is

[TABLE]

where the meaning of the integral as ordinary integral is justified because the integrand is a function. Rewriting the integrand, one obtains

[TABLE]

On the other hand, for all $x\not\in\{Z\}$,

[TABLE]

where the integral of the tempered distribution $f_{Z}$ is deemed defined as in appendix A. For controlling the difference it therefore suffices to bound the terms

[TABLE]

for various $z\in\{Z\}$. To this end, the elementary integral

[TABLE]

is further examined.

B.2 Analysis of the integral of $R_{N,T}$

It exhibits $R_{N,T}$ an oscillatory behaviour with an amplitude envelope which does not recede to zero (as $N\rightarrow\infty$), no matter at which place $x$ this property is considered. (This is in contrast to the $\text{si}(\cdot)$ function.) However, since $R_{N,T}$ bounded by this envelope $T^{-1}\cdot(\sin(2\pi\xi/(2T)))^{-1}$ for all $\xi\not\in T\cdot\mathbb{Z}$, its integral from [math] to $x$ tends pointwise (for these $x$) to a function, below seen to be $\tfrac{1}{2}+\left\lfloor x/T\right\rfloor$. The following shows that in fact the convergence is uniform on any closed interval within $[-T/2,T/2]$ not containing the zero. The main tool employed for this is standard Fourier analysis.

Let $J_{N,T}$ be as above. Since $R_{N,T}$ is alternatively written as $R_{N,T}(x)=T^{-1}\cdot(1+2\sum_{k=1}^{N-1}\cos(2\pi k\cdot x/T))$, it is $J_{N,T}$ equal to

[TABLE]

In order to show $J_{N,T}(x)\rightarrow\tfrac{1}{2}+\left\lfloor x/T\right\rfloor$ pointwise, define the periodic odd function

[TABLE]

with $h_{T}(x)$ (arbitrarily) set equal to zero otherwise. It suffices then to show that

[TABLE]

For this, one recognizes that $A_{N,T}$ is the truncated Fourier series expansion of $h_{T}$. (Set $T=1$ and employ the orthonormal system $\{v_{k}\}$ with $v_{k}(x^{\prime}):=\sqrt{2}\cdot\sin(2\pi k\cdot x^{\prime})$ on the interval $[0,1]$.) Standard Fourier analytical techniques then immediately yield convergence in the $L_{2}$-norm on $[0,T]$, since $h_{T}$ is piecewise continuous.

For obtaining the uniform convergence and a non-asymptotic (numerical) bound, a result building on the special form of the summands in $A_{N,T}$ can be employed. (See for example Theorem 6.5 (Abel’s test) in [Wal88], repeated below for convenience.)

Theorem 3:

(condensed from [Wal88]) Let $(f_{n})_{n\in\mathbb{N}}$ be a sequence of real-valued functions defined on $[a,b]$, fulfilling that the absolute value of the partial sums of $\sum_{n=1}^{\infty}f_{n}(x)$ are uniformly bounded, say by constant $M\in\mathbb{R}$. Let $(a_{n})_{n\in\mathbb{N}}\subset\mathbb{R}$ be a decreasing sequence converging to zero. Then $\sum_{n=1}^{N}a_{n}f_{n}(x)$ converges uniformly on $[a,b]$ as $N\rightarrow\infty$, and the residual sums fulfill for all $N\in\mathbb{N}$

[TABLE]

uniformly in $x$. (The proof of the theorem is based on Abel’s lemma.) ∎

The application of the theorem to the present case allows to derive a bound on the error caused by truncation of the infinite series, yielding the following result.

Lemma 3:

Let $A_{N,T}$ be defined as above. Let $I_{2}$ be a closed interval $\subset\mathbb{R}$ with $\{\left\lvert x-jT\right\rvert,x\in I_{2},j\in\mathbb{Z}\}\geq T\cdot\epsilon^{*}>0$, i.e. $I_{2}$ is bounded away from $T\cdot\mathbb{Z}$ by at least $T\cdot\epsilon^{*}$. Then $A_{N,T}(x)\rightarrow h_{T}(x)$ uniformly on $I_{2}$, and for all $x\in I_{2}$

[TABLE]

Proof: Similarly as in the proof of Theorem 4, an expression for $Q_{1}:=\sum_{k=1}^{N-1}\sin(2\pi k\cdot x/T)$ can be obtained. With $x^{\prime}=x/T$ and $0<x^{\prime}<1/2$

[TABLE]

With $\sin(u)\geq\tfrac{2}{\pi}u$ for $u\in[0,\pi/2]$ get

[TABLE]

(This bound is extended to $x^{\prime}\in(0,1)$, by observing the symmetry of $\sin(2\pi\cdot x^{\prime}/2)$, yielding $Q_{1}\leq 1/(2\min(x,1-x))$, and analogously to $x\in\mathbb{R}\backslash(T\cdot\mathbb{Z})$.) For $x\in I_{2}$, it follows

[TABLE]

Applying theorem B.2 yields the result.∎

Remark 3:

Since the above arguments found that

[TABLE]

the previous result becomes applicable on

[TABLE]

The left-hand side in turn equals $\left\lvert J_{N,T}(x)-(-\frac{1}{2}+1_{0\leq x/T})\right\rvert$ for $x\in(-T,T)$. ∎

In the above, the factor $T$ appearing next to $\epsilon^{*}$ reflects the fact that as we arbitrarily scale $Z$ (and let $\kappa$ be constant), the minimum distance of $I_{2}$ which needs to be kept from $T\cdot\mathbb{Z}$ must vary linearly in $T$ in order to leave the bound for $\left\lvert Q_{1}\right\rvert$ invariant. When employing the lemma, $\epsilon^{*}$ will be chosen (as usual) in dependence of $N$ such that the contribution to the total error from $J_{N,T}$ which turn out to be evaluated in the $[-\epsilon^{*},\epsilon^{*}]$ interval just balances with the contribution from the bound (55). I.e., $\epsilon^{*}$ may not be choosen too small.

In the following theorem, the ”max” expression appearing equals $\max_{z_{0}}P^{Z}([z_{0}-T\epsilon^{*},z_{0}+T\epsilon^{*}])$, where $P^{Z}$ is the measure of $Z$ derived from $\hat{F}_{n}$.

Theorem 4:

Let $X_{i},Z,F_{Z},T_{Z}$ be as in section 2 and $\tilde{F}$ be as in equation (21). Let $T>T_{Z}$ and $\kappa=T/T_{Z}$, and $x_{0}$ chosen with $\kappa\cdot z_{min}\leq x_{0}<z_{min}$. Let $M_{2}\in\mathbb{R}$ with

[TABLE]

Then

[TABLE]

for each $x\in I\backslash A_{d}$, where $A_{d}$ is the set of places of discontinuity of $F_{Z}$.

Proof: Let $x\in I\backslash A_{d}$. It is, with $f_{Z}=\tfrac{1}{N_{Z}}\sum_{z\in\{Z\}}\delta_{z}$,

[TABLE]

Since $x_{0}<\min Z$ and $\max Z-x_{0}<\kappa\cdot z_{max}-x_{0}\leq T$, it is $x_{0}-z<0$ and $x_{0}-z>x_{0}-\max Z>-T$, i.e. $x_{0}-z\in(-T,0)$ for all $z\in\{Z\}$. Therefore the lower evaluation of $J_{N,T}$ tends to $-\tfrac{1}{2}$ as $N\rightarrow\infty$. With $x\not\in A_{d}$, it is $x\neq z$. Via $x\in I$, it is $x-z\in(-T,T)$. For $x>z$, it tends $J_{N,T}(x-z)$ towards $\tfrac{1}{2}$. With this

[TABLE]

where the remark B.2 was used observing $x-z\in(-T,T)$.

Examine the second term first. For summands with $z$ close to $x$, Lemma B.2 will not be applicable. Thus let $\epsilon^{*}>0$ and represent $\{Z\}$ as union of the multisets $A_{1}(x):=\{z\in\{Z\},\left\lvert x-z\right\rvert\leq T\cdot\epsilon^{*}\}$ and $A_{2}(x):=\{z\in\{Z\},\left\lvert x-z\right\rvert>T\cdot\epsilon^{*}\}$. Then using the assumption (58) in the summation over $A_{1}$ yields

[TABLE]

Herein, it was used that $\left\lvert J_{N,T}(x-z)+\tfrac{1}{2}-1_{z\leq x}\right\rvert$ is bounded by $\tfrac{1}{2}$ uniformly in $N$ and $x$. Similarly, the first term in eqn. (65) is bounded, yielding the same bound. 666In fact, if $\kappa$ is sufficiently large or $\epsilon^{*}$ sufficiently small such that $T\cdot\epsilon^{*}<(\kappa-1)z_{min}$, then $A_{1}(x_{0})$ will be empty, which effects below a reduction in the ultimate bound by factor $1/\sqrt{2}$. In total

[TABLE]

It is easy to determine that the optimum choice for $\epsilon^{*}$ minimizing $E$ yields

[TABLE]

which proves the theorem. ∎

The overall result obtained therefore is: For places $x$ in the discontinuity set $A_{d}$ of $F_{Z}$, no statement is made. For places $x$ outside this set, Theorem B.2 states a bound on the pointwise error which uniformly holds over $I\backslash A_{d}$.

Plausible values for $M_{2}$ can be only estimated crudely before having obtained an estimate of $f_{Z}$. The following remark provides a hint to possible occuring values in the non-degenerate case (usually resulting from $X$ having a continuous distribution) by deriving an $M_{2}$ for the normal distribution. In cases where the sample $(X_{i})_{i=1...n}$ exhibits regularities (e.g. $X_{i}\in\mathbb{Z}$ for all $i$) or even is degenerate ($X_{i}=\text{const}$), the $M_{2}$ will have excessive values.

Remark 4: Let $U\sim\mathcal{N}(0,\sigma^{2})$ and $T=5\sigma$. Let $0<\epsilon^{*}<5$. Then

[TABLE]

i.e. for the random variable $U$ the appropriate value for $M_{2}$ would be $10/\sqrt{2\pi}$. ∎