hermiter: R package for Sequential Nonparametric Estimation

There are existing algorithms for online estimation of all of the aforementioned statistical quantities. For PDFs, existing algorithms include recursive kernel estimators as discussed in Chapters 4 and 5 of Greblicki and Pawlak (2008) and Chapter 7 of Devroye and Gyorfi (1985). See also Slaoui and Jmaei (2019) for recursive estimators based on Bernstein polynomials. In the case of CDFs, existing algorithms include those based on the empirical distribution function estimator, the smooth kernel distribution function estimator (Slaoui, 2014) and Bernstein polynomials (Jmaei et al, 2017). These PDF and CDF estimators furnish sequential estimates of the probability density and cumulative probability at predefined values of the support of the probability density function, which is not as general as the online estimation of the full PDF and CDF.

Several algorithms exist for estimating quantiles in a sequential (one-pass) manner. These algorithms can be differentiated into two broad categories. The first category is comprised of those algorithms that directly approximate arbitrary empirical quantiles of a set of numeric values using sub-linear memory. Here, the empirical p-quantile is defined as $x_{(\lceil pn\rceil)}$ for $n$ observations, where $x_{(i)},\,i\in 1,\dots n$, are the observations sorted in ascending order. No particular data generating process is presupposed. These algorithms provide either deterministic or probabilistic guarantees on the rank error of approximate sample quantiles and have space requirements that grow with required accuracy and potentially also with the number of observations. Noteworthy algorithms in this category include those due to Greenwald and Khanna (2001) and Karnin et al (2016). For recent reviews of algorithms of this type see Greenwald and Khanna (2016), Luo et al (2016) and Chen and Zhang (2020). Rank error accuracy is but one choice of quality metric however and it has been noted that controlling rank error accuracy does not preclude arbitrarily large relative error (for heavily skewed, or heavy tailed data) as per Masson et al (2019). In Masson et al (2019) an algorithm providing relative error accuracy guarantees for estimating sample quantiles (with space requirements that grow with accuracy) is proposed to address this (see also Epicoco et al (2020) for an enhancement of this algorithm). It is also worth noting that many of the aforementioned algorithms cannot be used to form dynamic quantile estimates based on a recent subset of observations since they are insertion-only algorithms. Here insertion-only means that the estimates can only be updated with new observations in an online manner, but previous observations cannot be removed in such a manner. Exceptions include those algorithms that also allow deletion operations and use turnstile semantics to maintain an estimate over a sliding window. These algorithms include the Dyadic Count-Sketch (Luo et al, 2016) for example. Such algorithms assume a fixed universe of potential values for the values being analyzed and typically have higher space and update time complexity than insertion-only algorithms. See also Zhao et al (2021) for recent developments using a bounded deletion model.

The second broad category of online quantile estimation algorithms, of which the Hermite series based algorithm is a member (i.e. the algorithm proposed in Stephanou et al (2017) and implemented in hermiter), assume that the observations to be analyzed are a realization of a random process. Approximating quantiles is then a nonparametric statistical estimation problem. Algorithms in this category are numerous (Jain and Chlamtac, 1985; Raatikainen, 1987, 1990; Naumov and Martikainen, 2007; Chen et al, 2000; Yazidi and Hammer, 2017; Hammer et al, 2019b, a; Tiwari and Pandey, 2019). The aforementioned algorithms in the literature only apply to pre-specified quantiles however. A recently introduced algorithm by Gan et al (2018) involving moment-based quantile sketches, furnishes an efficient approach to quantile estimation but focusses on rapid mergeability and single quantile queries. See also Mitchell et al (2021) for a comparison of the moment-based quantile sketch algorithm with several other algorithms on a large collection of real-world data sets. In Mitchell et al (2021) orthogonal series based quantile estimators similar to those in Stephanou et al (2017) are studied, including an estimator based on the Chebyshev-Hermite polynomials (also known as the “Probabilist’s" Hermite polynomials). These are different orthogonal polynomials than those used in Stephanou et al (2017), which are instead the “Physicist’s" Hermite polynomials, and thus the results in Mitchell et al (2021) do not directly apply to the estimators implemented in hermiter. The popular t-digest approach (Dunning, 2021) maintains an estimate of the empirical distribution function and allows the estimation of arbitrary quantiles, presenting similar capabilities to the Hermite series based approach. This algorithm has the useful property of maintaining an accuracy which is nearly constant relative to $p(1-p)$ for the p-quantile. In addition, the t-digest algorithm has many desirable properties in practice, such as approximate mergeability, and is fast and accurate.

In the bivariate setting, there is only one alternate approach to online estimation of nonparametric correlation coefficients known to the authors, namely the algorithms introduced in Xiao (2019) which involve discretizing the joint distribution of the random variables.

The Hermite series density estimator for the univariate PDF, $f(x)$, is defined as follows (Schwartz, 1967; Walter, 1977; Greblicki and Pawlak, 1984, 1985; Liebscher, 1990),

where $\mathbf{x_{i}}\sim f(x),\,i=1,2,\dots n$, and $h_{k}(x)=\left(2^{k}k!\sqrt{\pi}\right)^{-\frac{1}{2}}e^{-\frac{x^{2}}{2}}H_{k}(x),\,k=0,1,\cdots,N$ are the normalized Hermite functions defined from the Hermite polynomials, $H_{k}(x)=(-1)^{k}e^{x^{2}}\frac{\mathrm{d}^{k}}{\mathrm{d}x^{k}}e^{-x^{2}},\,k=0,1,\cdots,N$.

The bivariate Hermite series density estimator for the bivariate PDF, $f(x,y)$, has the following form:

where $(\mathbf{x_{i}},\mathbf{y_{i}})\sim f(x,y),\,i=1,2,\dots n$. Note that although standardization of the observations is not required for Hermite series based estimators, our empirical studies have revealed that it often improves accuracy. Here standardization refers to subtracting the estimated mean from each observation and dividing by the estimated standard deviation. The estimates for the Hermite series coefficients based on these standardized observations are biased compared to the coefficients that would have been obtained by using the true mean and standard deviation for standardization. Despite this bias, we have typically observed better overall results in practice using standardization. The PDF estimators previously defined need to be modified for standardized observations. For the univariate case,

where $\hat{\sigma}$ is the estimated standard deviation. For the bivariate case,

where $\hat{\sigma}^{(1)}$ and $\hat{\sigma}^{(2)}$ are the estimates of the standard deviation of $x$ and $y$ respectively.

The univariate density estimator (1) can be used to define the following univariate CDF estimator,

Quantiles can be obtained numerically by applying a root-finding algorithm to solve for the quantile value, $\hat{x}_{p}$, at cumulative probability $p$ via,

In Stephanou et al (2017) an alternate univariate CDF estimator was also introduced,

We have found that while the asymptotic properties of the alternate univariate CDF estimator (4) are the same as (3), the empirical, finite sample performance of (4) is superior for quantile estimation specifically, see Stephanou et al (2017). A vectorized implementation of the bisection root finding algorithm is a robust way to calculate the quantiles at an arbitrary vector of cumulative probability values. This is one of the two algorithms available for quantile estimation in hermiter. The second algorithm allows for more rapid but sometimes less accurate approximation of the quantiles by calculating the CDF values at a fixed vector of values and applying linear interpolation to find the vector of quantiles using the R function findInterval. It is also noteworthy that in (3) and (4), it is computationally advantageous to utilize the recurrence relations for integrals of Hermite functions in (5) and (6),

and

where $\mbox{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-t^{2}}dt$ is the complementary error function. Convergence acceleration techniques can be successfully applied to improve accuracy at a given order of truncation for Hermite series. In particular, hermiter applies straightforward iterated averaging techniques to the truncated series (see Boyd and Moore (1986) and Table A.2 of Boyd (2018)) to improve the accuracy of univariate PDF and CDF estimates along with quantile estimates obtained through the CDF estimates. We have observed that empirically, applying series acceleration typically yields better results. Deriving formal guarantees that these methods are generally effective in our context is an area for future research. It is worth noting that the results presented in Stephanou et al (2017) can be reproduced by setting the algorithm used for quantile estimation to ‘bisection’ and the accelerate series argument to FALSE.

In a similar manner to the univariate case, the bivariate density estimator (2) can be used to construct the following bivariate CDF estimator:

The recurrence relation in (5) can again be applied for rapid calculation of the integrals.

We can also use the Hermite series based estimators defined above to obtain estimators for the Spearman Rho and Kendall Tau nonparametric correlation measures. For the Spearman correlation coefficient estimator, $\hat{R}_{N}$, we plug the bivariate Hermite PDF estimator and univariate Hermite CDF estimator into a large sample form of the Spearman correlation coefficient (see Stephanou and Varughese (2021b)) to obtain,

where $\hat{F}_{N,n}^{x}(x),\hat{F}_{N,n}^{y}(y)$ correspond to univariate CDF estimators as defined in equation (3), updated on $x$ and $y$ observations respectively. The function $\hat{f}_{N,N,n}(x,y)$ is the bivariate PDF estimator defined in equation (2). We can phrase the estimator above in terms of linear algebra operations for computational efficiency:

where $\mathbf{\hat{a}_{(1)}}^{(n)}$ are the coefficients $(\mathbf{\hat{a}_{(1)}}^{(n)})_{k}$ associated with the Hermite CDF estimator $\hat{F}_{N,n}^{x}(x)$, $\mathbf{\hat{a}_{(2)}}^{(n)}$ are the coefficients $(\mathbf{\hat{a}_{(2)}}^{(n)})_{k}$ associated with the Hermite CDF estimator $\hat{F}_{N,n}^{y}(y)$ and $\mathbf{\hat{A}}^{(n)}$ are the coefficients $(\mathbf{\hat{A}}^{(n)})_{kl}$ associated with the Hermite bivariate PDF estimator $\hat{f}_{N,N,n}(x,y)$. The matrix $(\mathbf{W})_{kl}=\int_{-\infty}^{\infty}h_{k}(u)\int_{-\infty}^{u}h_{l}(v)dvdu$ and vector $(\mathbf{z})_{k}=\int_{-\infty}^{\infty}h_{k}(u)du$. In all the matrices and vectors we have defined, $k,l=0,\dots,N$. These integrals can be evaluated numerically once and the values stored for rapid calculation. The sequential nature of the estimators follows from the fact that we have obtained sequential update rules for the coefficients and the fact that there is no explicit dependence on the observations in these expressions as discussed in section 6.4 below.

For the Kendall Tau estimator, $\hat{\tau}_{N}$, we plug the bivariate Hermite PDF estimator and bivariate Hermite CDF estimator into a form of the Kendall correlation coefficient relating directly to the probability of concordance (see Gibbons and Chakraborti (2010) for more details on this representation), to obtain,

where $\hat{F}_{N,N,n}(x,y)$ is the bivariate CDF estimator defined in equation (7) and $\hat{f}_{N,N,n}$ is the bivariate PDF estimator defined in equation (2). We can again phrase the estimator above in terms of linear algebra operations for computational efficiency:

The particular utility of all the Hermite series based estimators in online/sequential estimation follows from the ability to update the estimates of the Hermite series coefficients in an online, $O(1)$ manner. In particular, for the coefficients in (1), we can define the following updating algorithm,

For the coefficients in (2), we can define an updating algorithm in terms of the outer product of $\mathbf{h}(\mathbf{x})$ and $\mathbf{h}(\mathbf{y})$,

In the update rules (10) and (11) above, $\mathbf{h}(\mathbf{x}_{i})$ and $\mathbf{h}(\mathbf{y}_{i})$ are the vectors $h_{k}(\mathbf{x}_{i}),\,k=0,\dots,N$ and $h_{l}(\mathbf{y}_{i}),l=0,\dots,N$ respectively. For computational efficiency it is advantageous to calculate $\mathbf{h}(\mathbf{x}_{i})$ and $\mathbf{h}(\mathbf{y}_{i})$ making use of the recurrence relation for the Hermite polynomials $H_{k+1}(x)=2xH_{k}(x)-2kH_{k-1}(x)$. The computational cost of updating the coefficients, $\mathbf{\hat{a}}^{(i)}$ and $\mathbf{\hat{A}}^{(i)}$ as above is manifestly constant ($O(1)$) with respect to the number of previous observations. In non-stationary scenarios the following exponential weighting scheme can be applied for the coefficients defined in (1),

where $0<\lambda\leq 1$ controls the weight of new terms associated with new observations. For the coefficients defined in (2) we have,

The computational cost of updating the coefficients, $\mathbf{\hat{a}}^{(i)}$ and $\mathbf{\hat{A}}^{(i)}$ is again manifestly constant ($O(1)$) with respect to the number of previous observations.

As noted above, the Hermite series based estimators have been found to typically perform better on standardized observations. The observations can be standardized in an online manner by subtracting an online estimate of the mean and dividing by an online estimate of the standard deviation.

We present a summary of the core functions provided by the hermiter package in Table 2 and elaborate on the functions in more detail in the remainder of the section.

A hermite_estimator R S3 object is constructed as below. The argument, N, adjusts the number of terms in the Hermite series based estimator and controls the trade-off between bias and variance. A lower N value implies a higher bias but lower variance and vice versa for higher values of N. We have found that empirically, different contexts benefit from different default values of N for good performance. We have distilled our experience on a wide array of data sets into defaults that have been implemented in the hermiter package. For univariate, non-exponentially weighted estimators, the default is N$=50$ whereas for univariate, exponentially weighted estimators, the default is N $=20$. For bivariate, non-exponentially weighted estimators, the default is N $=30$. Finally, for bivariate, exponentially weighted estimators, the default is N $=20$. For further detail on selecting $N$ in the univariate setting, see section 4.1 of Stephanou et al (2017) and for the bivariate setting, see section 3.3 of Stephanou and Varughese (2021b). The argument, standardize, controls whether or not to standardize observations before applying the estimator. Standardization usually yields better results and is recommended for most estimation settings, thus the default for standardize is set to TRUE. Finally, the estimator can be optionally initialized with a batch of observations via the argument observations.

A univariate estimator is constructed as follows (the default estimator type is univariate, so this argument does not need to be explicitly set):

Similarly for constructing a bivariate estimator:

If one wishes to initialize the estimator with an initial batch of observations, the following syntax can be utilized. For univariate observations we have:

For bivariate observations we have:

Note that the constructors above only take numeric observations arguments i.e. a numeric vector or numeric matrix respectively. Initializing with data frames is not presently supported even if all variables are numeric.

In the sequential setting, observations are revealed one at a time. A hermite_estimator object can be updated sequentially with a single or multiple new observations by utilizing the update_sequential method. Note that when updating the Hermite series based estimator sequentially, observations are also standardized sequentially if the standardize argument is set to TRUE in the constructor. This means that using the update_sequential method with standardize set to TRUE for an initial set of observations will give a slightly different result to initializing the estimator using the observations argument. For updating with univariate observations we have:

For bivariate observations we have:

One can verify that the new observations have been incorporated by printing the hermite_estimator object before and after the sequential updates: the number of observations can be seen to have increased by the appropriate amount. The generic method update_sequential dispatches to the appropriate implementation depending on whether the hermite_estimator object is univariate or bivariate. In the univariate case, the sequential updating of the Hermite series coefficients is based on (10). In the bivariate case, the sequential updating of the Hermite series coefficients is based on (11).

The hermiter package is also applicable to non-stationary data streams. An exponentially weighted form of the Hermite series coefficients can be applied to handle this case. The estimators will adapt to the new distribution and incrementally “forget" the old distribution. In order to use the exponentially weighted form of the hermite_estimator, the exp_weight_lambda argument must be set to a non-NA value. Typical values for this parameter are 0.01, 0.05 and 0.1. The lower the exponential weighting parameter, the slower the estimator adapts and vice versa for higher values of the parameter. However, the variance of estimands obtained from the hermite_estimator object increases with higher values of exp_weight_lambda, so there is a trade-off to bear in mind. For univariate observations we have:

For bivariate observations we have:

Again the generic method update_sequential dispatches to the appropriate implementation depending on whether the hermite_estimator object is univariate or bivariate. In the univariate case, the sequential updating of the Hermite series coefficients is based on (12). In the bivariate case, the sequential updating of the Hermite series coefficients is based on (13).

Hermite series based estimators can be consistently merged in both the univariate and bivariate settings as discussed previously. In particular, when standardize = FALSE, the results obtained from merging distinct hermite_estimator objects updated on subsets of a data set are exactly equal to those obtained by constructing a single hermite_estimator and updating on the full data set (corresponding to the concatenation of the aforementioned subsets). This holds true for the PDF, CDF and quantile results in the univariate case and the PDF, CDF and nonparametric correlation results in the bivariate case. When standardize = TRUE, the equivalence is no longer exact, but is accurate enough to be practically useful. Merging hermite_estimator objects is illustrated below. Note that the merge_hermite function takes in a list of hermite_estimator objects as an argument. In order to be merged, these objects must have a consistent type, either "univariate" or "bivariate". As an example in the univariate case:

The above example yields the following output which verifies that the merged hermite_estimator i.e. hermite_est_merged is exactly equal equal to the full hermite_estimator, i.e. hermite_est_full when standardize = FALSE.

A univariate example with standardize = TRUE is the following:

This yields the output:

The above example illustrates that the merged hermite_estimator i.e. hermite_est_merged is a close approximation to the full hermite_estimator, i.e. hermite_est_full when standardize = TRUE but it is not exactly the same object. As an example for the bivariate case:

When standardize = FALSE, merging utilizes (14) and (15) to combine the univariate and bivariate Hermite series coefficients respectively. When standardize = TRUE, merging utilizes the methods in Schubert and Gertz (2018) to combine the mean and standard deviation estimates of the individual hermite_estimator objects along with (16) and (17) to combine the univariate and bivariate Hermite series coefficients respectively.

A few final notes are that the standardize parameter must be the same i.e. either all TRUE or all FALSE for the list of hermite_estimator objects to be merged. In addition, the N argument must be the same for all estimators in order to merge them.

The main advantage of Hermite series based estimators is that they can be updated in a sequential/one-pass manner as above and subsequently probability densities and cumulative probabilities at arbitrary x values can be obtained, along with arbitrary quantiles. The hermite_estimator object only maintains a small and fixed number of coefficients and thus uses minimal memory. The syntax to calculate probability densities, cumulative probabilities and quantiles in the univariate setting is presented below. The PDF estimators (1) and (2) are implemented in the dens method, where dispatch occurs to the appropriate implementation depending on whether the hermite_estimator object is univariate or bivariate. As described in section 6.2.1, density estimates can be negative in principle. The parameter clipped can be set to TRUE in order to replace negative values with a small, but strictly positive value (i.e. 1e-8).

The CDF estimators (3) and (7) are implemented in the cum_prob method, where again dispatch occurs to the appropriate implementation depending on whether the hermite_estimator object is univariate or bivariate. Cumulative probability estimates may not be monotonically non-decreasing and might lie outside the range $[0,1]$ as per section 6.2.1. These limitations can be partially addressed by clipping the estimates to lie between $[0,1]$, this can be done by setting the clipped parameter to TRUE. There are also various post-processing strategies that can be employed to enforce the monotone property (e.g. making use of the R function cummax). The quantile function quant makes use of the CDF estimator (4) and is only defined for univariate hermite_estimator objects. It is worth noting that since the Hermite series based quantile estimator is directly estimating the quantile function and not the empirical quantiles, we do not expect the 0-quantile to correspond to the minimum observation, nor do we expect the 1-quantile to correspond to the maximum observation. In fact, since the Hermite series based estimators are defined on the full real line, the 0-quantile (true minimum) is $-\infty$ and the 1-quantile (true maximum) is $\infty$ by definition. As a result, estimating these quantiles, i.e. estimating the minimum and maximum directly is not sensible using these estimators. Quantiles for $p$ very close to 0 and $p$ very close to 1 may still be reasonably accurate however. The default behaviour of the quant function is therefore to return a quantile for $p$ very close to 0 for the 0-quantile and a quantile for $p$ very close to 1 for the 1-quantile.

As an example of estimating the PDF, CDF and quantiles in the univariate case:

Actual values of the PDF, CDF and quantiles are evaluated as below:

The estimated versus actual values are visualized in Figures 1, 2 and 3.

The aforementioned suitability of Hermite series based estimators in sequential and one-pass batch estimation settings extends to the bivariate case. Probability densities and cumulative probabilities can be obtained at arbitrary points. The syntax to calculate probability densities and cumulative probabilities along with the Spearman and Kendall correlation coefficients in the bivariate setting is presented below. Note that the functions spearmans and kendall implement the equations (8) and (9) respectively. As an example of estimating the PDF, CDF and nonparametric correlation coefficients in the bivariate case:

Actual values of the bivariate PDF, CDF, Spearman and Kendall correlation coefficients are evaluated as below.

The estimated versus actual values for the PDF and CDF are visualized in Figures 4, 5 respectively.

The Spearman and Kendall correlation coefficient results are presented in Table 3.

In the univariate case, simulation studies have been conducted to demonstrate competitive performance of the Hermite series based algorithms in the sequential quantile estimation setting in Stephanou et al (2017). In addition, simulation studies have been performed to evaluate the accuracy of the Hermite series based CDF estimator in Stephanou and Varughese (2021a) where it was concluded that while smooth kernel distribution function estimators are more accurate in the non-sequential setting, the performance of the sequential Hermite series based CDF estimators is still sufficiently good to be practically useful in online settings. The recent tdigest package in R offers very similar functionality for sequential quantile estimation to hermiter. This similarity, along with the availability of a R implementation and the popularity of the tdigest algorithm suggests a new and useful comparison to the Hermite series based algorithms in terms of accuracy and performance in the stationary setting. Note that in the non-stationary setting, where quantiles may vary over time, only the hermiter package can be directly applied, implying that even if the accuracy and performance of tdigest and hermiter are similar in the stationary setting, the hermiter package still adds valuable functionality. The simulation comparison we conduct utilizes the suite of test distributions provided by the R package benchden (Mildenberger and Weinert, 2012) excluding heavy-tailed distributions. The reason we exclude the heavy-tailed densities from the benchden test suite is that most of these have undefined variances and some even have undefined means. The Hermite series based quantile estimators as implemented in hermiter require at least the first two moments to exist and thus we cannot meaningfully compare hermiter and tdigest on these densities. The simulation study is as follows: for sample sizes, $n=10,000\mathchar 24635\relax\;100,000\mathchar 24635\relax\;1000,000\mathchar 24635\relax\;10,000,000$ and each of the test distributions, the following steps are repeated.

1. 1.

   For each test distribution, we initialize a hermite_estimator object (using defaults N $=50$, standardize = TRUE), with $n$ observations drawn from the test distribution. We record the integrated absolute error (IAE) and partial integrated absolute error (pIAE) between estimated quantiles - obtained via the quant function with algorithm = "interpolate" and accelerate_series = TRUE - and the true quantile values. Similarly for the tdigest algorithm using the tdigest and quantile functions. The IAE is defined as,

   $$\mbox{IAE}=\int_{0}^{1}\lvert\hat{q}(p)-q(p)\rvert dp,$$

   and the partial IAE is defined as,

   $$\mbox{pIAE}=\int_{0.01}^{0.99}\lvert\hat{q}(p)-q(p)\rvert dp,$$

   where $\hat{q}(p)$ is the estimated quantile value at cumulative probability $p$ and $q(p)$ is the true quantile value at cumulative probability $p$. The IAE captures the integrated absolute error across the full range of cumulative probability values. The partial IAE captures the integrated absolute error excluding extreme quantiles (namely those where $p<0.01$ and $p>0.99$). In practice, we apply Quasi-Monte Carlo integration with a Sobol sequence to numerically determine the value of these integrals.

2. 2.

   We repeat the above $m=100$ times and average across all runs to obtain estimates of the mean integrated absolute error, $\mbox{MIAE}=E\left(\mbox{IAE}\right)$, and partial mean integrated absolute error, $\mbox{pMIAE}=E\left(\mbox{pIAE}\right)$.

This yields a MIAE and pMIAE value for each test distribution and sample size $n$ for both the Hermite series based quantile estimator and tdigest based quantile estimation algorithm. There are $21$ non heavy-tailed test distributions in total. We summarize the results of comparing the MIAE of hermiter and tdigest estimates in Table 5. The results for pMIAE are presented in Table 6.

The results reveal that the accuracy of the Hermite series based approach is better than the tdigest algorithm for fairly large numbers of observations even for the full MIAE error measure. In this simulation study we are particularly interested in performance on a number of common distributions across a range of non-extreme (i.e. more central) quantiles, for which the pMIAE results are more relevant. These results demonstrate significantly better accuracy using the Hermite series based algorithm for quantile estimates formed from more than $n=10,000$ observations. The distributions in benchden include those with full real-line, half real-line and compact support and thus the simulation study should yield information about typical performance in a broad range of scenarios. We reiterate that we exclude the heavy-tailed densities from the benchden test suite but should mention that the tdigest package appears to perform well in quantile estimation on all these heavy-tailed densities. Thus, for applications involving heavy-tailed distributions, tdigest is likely better suited.

In terms of computational performance, we plot the time taken to initialize the hermiter estimator and the tdigest estimator with $n=1000,000$ observations in Figure 8. It is apparent that hermiter is faster across all parameters of N considered. This is mainly the result of leveraging multithreaded computation on multiple cores for batch updating. This is enabled by default. Multithreaded computation can be disabled by setting options(hermiter.parallel = FALSE). Computation time will typically increase on multi-core systems in this case however, compare Figure 8 (multithreaded computation enabled) to Figure 9 (multithreaded computation disabled) for example. Both benchmarks were conducted on the same eight core machine (see section Computational details).

Finally, we plot the time taken to calculate vectors of quantile values ranging in length from a single quantile to 100,000 quantiles in Figure 10. The performance of tdigest is better for smaller numbers of quantiles but the distributions of run-times of hermiter and tdigest are similar for larger numbers of quantiles. Calculating large vectors of quantiles simultaneously is meaningful for drawing random variables for example. Given the fact that further implementation optimization is possible for hermiter (and potentially tdigest) and that performance results are likely operating system dependent, we regard the computational performance as roughly similar (while acknowledging that tdigest appears faster for smaller vectors of quantiles, perhaps due to less computational overhead).

In the bivariate case, extensive simulation studies have been conducted to demonstrate competitive performance in the sequential estimation of the Spearman correlation coefficient versus the only known competing algorithm in Stephanou and Varughese (2021b). In this section we present a succinct simulation study using essentially the same methodology as Stephanou and Varughese (2021b) where we extend the results to the sequential estimation of the Kendall Tau coefficient. We term the competing algorithms in Xiao (2019) count matrix based algorithms.

The simulation study involves generating i.i.d. samples of increasing size, $n=10,000\mathchar 24635\relax\;50,000\mathchar 24635\relax\;100,000$, from a bivariate normal distribution with mean vector, $\mu=(0,0)$, and covariance matrix, $\Sigma=(\sigma_{1},\rho\sigma_{1}\sigma_{2}\mathchar 24635\relax\;\rho\sigma_{1}\sigma_{2},\sigma_{2})$, where $\sigma_{1}=1,\sigma_{2}=1$. The correlation of the bivariate normal distribution is varied by drawing samples across a range of correlation parameters, $\rho=-0.75,-0.5,-0.25,0.25,0.5,0.75$. For each value of the sample size $n$ and correlation parameter, $\rho$, the below steps are repeated $m=100$ times.

1. 1.

   Draw a sample of n i.i.d. observations, $(\mathbf{x_{i}},\mathbf{y_{i}}),\quad i=1,\dots,n$ from the bivariate normal distribution with correlation parameter $\rho$, mean vector $\mu$, and covariance matrix $\Sigma$, using the rmvnorm function in the package mvtnorm (Genz et al, 2020; Genz and Bretz, 2009).

2. 2.

   Construct a bivariate hermite_estimator using hermiter and initialize with the full sample.

3. 3.

   Estimate the Spearman Rho and Kendall Tau correlation coefficients for the
   hermite_estimator object using the spearmans and kendall functions in hermiter.

4. 4.

   Estimate the true Spearman Rho on the full sample using the R function cor. Calculate the true Kendall Tau analytically using the relation $\tau=\frac{2}{\pi}\arcsin\rho$ which applies to bivariate normal distributions (Gibbons and Chakraborti, 2010). We use the exact analytical Kendall Tau since the implementation in cor is prohibitively slow for estimating Kendall Tau with larger sample sizes (as it would be in other standard implementations).

5. 5.

   Calculate the absolute error between the Hermite series based Spearman Rho and Kendall Tau estimates and the respective exact values.

The mean absolute error (MAE) between the Hermite series based estimates and the exact Spearman and Kendall rank correlation coefficients for a particular value of $\rho$ and sample size $n$ are then estimated through $\mbox{MAE}(\hat{R}_{N})=\frac{1}{m}\sum_{j=1}^{m}\lvert\hat{R}^{(j)}_{N}-R^{(j)}\rvert$ and $\mbox{MAE}(\hat{\tau}_{N})=\frac{1}{m}\sum_{j=1}^{m}|\hat{\tau}^{(j)}_{N}-\tau^{(j)}\rvert$ respectively, where $j$ indexes a particular set of $n$ observations. The same steps as above are repeated for the count matrix algorithms namely, algorithms 2 and 3 of Xiao (2019). For the choice of parameters, we use the default of hermiter for non-exponentially weighted, bivariate estimators, N $=30$, which we have found to give good performance in a wide range of scenarios. For algorithm 2 of Xiao (2019) we use the recommended value of 30 cutpoints for estimating Spearman Rho. For algorithm 3 of Xiao (2019) we use the recommended value of 100 cutpoints for estimating Kendall Tau. It is worth noting that the number of values to be maintained in memory is significantly larger in the count matrix case with 100 cutpoints than the Hermite estimator case with N $=30$. The results are summarized as the average MAE (and standard deviation of MAE) across all the values of $\rho$ in Table 7 for the Hermite series based algorithm with N $=30$ and count matrix algorithm with $c=30$ for Spearman Rho estimation. Similarly, results are summarized as the average MAE (and standard deviation of MAE) across all the values of $\rho$ in Table 8 for the Hermite series based algorithm with N $=30$ and count matrix algorithm with $c=100$ for Kendall Tau estimation. It is clear that the Hermite series based algorithms present an accuracy advantage. The mean MAE and variation in MAE are also summarized in Figure 11.

The computational performance of the algorithms in hermiter is roughly comparable with the algorithms in Xiao (2019). That said, when updating the Spearman or Kendall correlation coefficient with each new observation, we have observed better computational performance with the Hermite series based algorithms. We acknowledge that this may depend on the implementation of the count matrix algorithms however. We omit detailed benchmarking for brevity. A final note is that the calculation time for the Hermite series based Kendall Tau estimator is thousands of times faster than the cor function in R with $n=10,000$ observations and this differential grows rapidly with $n$. For large samples (e.g. millions or even billions of observations), the standard methods of calculating the Kendall correlation coefficient are intractable. The Hermite series based approach is not only tractable but very fast in practice.