A comparison of efficient approximations for a weighted sum of chi-squared random variables

Inmany applications, the cumulative distribution function (cdf) FQN of a positively weighted sum of N i.i.d. chi-squared randomvariables QN is required.Although there is no known closed-form solution for FQN , there are many good approximations. When computational efficiency is not an issue, Imhof’smethod provides a good solution. However, when both the accuracy of the approximation and the speed of its computation are a concern, there is no clear preferred choice. Previous comparisons between approximatemethods could be considered insufficient. Furthermore, in streaming data applications where the computation needs to be both sequential and efficient, only a few of the available methods may be suitable. Streaming data problems are becoming ubiquitous and provide the motivation for this paper. We develop a framework to enable a much more extensive comparison between approximate methods for computing the cdf of weighted sums of an arbitrary random variable. Utilising this framework, a new and comprehensive analysis of four efficient approximate methods for computing FQN is performed. This analysis procedure is much more thorough and statistically valid than previous approaches described in the literature. A surprising result of this analysis is that the accuracy of these approximate methods increases with N . Electronic supplementary material The online version of this article (doi:10.1007/s11222-015-9583-4) contains supplementary material, which is available to authorized users. B Dean A. Bodenham d.bodenham10@imperial.ac.uk 1 Department of Mathematics, Imperial College London, London, UK 2 D-BSSE, ETH Zürich, Basel, Switzerland 3 Heilbronn Institute of Mathematics, University of Bristol, Bristol, UK


Introduction
The cumulative distribution function (cdf) F Q N of a positively weighted sum of i.i.d. χ 2 1 random variables Q N , has no known closed-form solution. An approximation of F Q N is used in goodness-of-fit tests (Moore and Spruill 1975) and various other applications (Zhang and Chen 2007;Jayasuriya 1996;Bentler and Xie 2000). Our particular interest is change detection in streaming data (Bodenham 2014).
In offline situations where computational resources are not an issue, Imhof's method (Imhof 1961), which inverts the characteristic function numerically, should be the preferred choice. It can be considered exact (Solomon and Stephens 1977;Johnson et al. 2002) since it provides error bounds and can be used to compute F Q N (x), for some quantile value x, to within a desired precision. Similar numerical methods such as Farebrother's method (Farebrother 1984) could also be used, but some (Sheil and O'Muircheartaigh 1977;Davis 1977;Davies 1980) lack the precision-bounding feature of Imhof's method. However, Imhof's method and Farebrother's method are both iterative, which affects their speed of computation, as shown in Sect. 6.4. Besides being iterative, these methods all require the entire vector of coefficients (d 1 , . . . , d N ) to be stored in order to compute the approximate cdf. As described in Sect. 2, this may not be possible in a streaming data context. Perhaps, the earliest approximate method, which has come to be known as the Satterthwaite-Welch method (Welch 1938;Satterthwaite 1946;Fairfield-Smith 1936), involved matching the first two moments of Q N with the first two moments of a Gamma distribution. See Box (1954, Sect. 3) for a discussion on the history of this method. The Hall-Buckley-Eagleson (Hall 1983;Buckley and Eagleson 1988) and Wood F (Wood 1989) methods match the first three moments of Q N to other distributions in a similar fashion. The Lindsay-Pilla-Basak method (Lindsay et al. 2000) matches the first 2n moments of Q N to a mixture distribution. These four moment-matching methods are described in Sect. 3 and are implemented in the R package momentchi2 (Bodenham 2015).
The method described in Solomon and Stephens (1977) takes the Satterthwaite-Welch method a step further by matching the first three moments of Q N to a random variable a X b , where X ∼ χ 2 1 . It is accurate in both the upper and lower tails, but requires the solution of two simultaneous non-linear equations, perhaps via an iterative method. An interesting method using Laguerre polynomials is described in Castaño-Martínez and López-Blázquez (2005), but is also iterative and requires the setting of certain control parameters.
While the methods discussed here have superseded those published previously (e.g. Patnaik 1949;Jensen and Solomon 1972), a good review of older methods can be found in Johnson et al. (2002). Although not considered here, a review of the current state-of-the art for weighted sums of non-central chi-squared random variables can be found in Duchesne and Lafaye De Micheaux (2010), and methods for computing the cdf of a single non-central chi-squared random variable are described in Farebrother (1987), Ding (1992) and Penev and Raykov (2000). An earlier version of this work appeared in the unpublished PhD thesis of Bodenham (2014).

Approximations in a streaming data context
If we wished to simply compute a single evaluation of F Q N , for some vector of coefficients d = (d 1 , d 2 , . . . , d N ), then we have already described a plethora of methods from which to choose. Amongst these, since Imhof's method is essentially exact it would probably be the preferred choice. There are, however, situations when Imhof's method might not be suitable. For instance, one might wish to compute F Q N (x), for Q N defined in Eq. (1), and then soon afterwards compute Imhof's method requires the whole vector of weights d in order to compute F Q N +1 (x ), but in a streaming data context (discussed in the next paragraph) N might be very large, and so storing the whole coefficient vector (d 1 , . . . , d N , d N +1 ) would be undesirable. Finally, Imhof's method is also iterative, since it runs until a specified precision is obtained. This is also unappealing, since iterative methods have the potential to be slow and computationally expensive. Given this construction, Imhof's method is clearly not suitable for deployment.
Streaming data algorithms (e.g. Gama et al. 2010;Bodenham and Adams 2013) require methods that are both fast and only require a small, fixed number of parameters and data to be stored. Amongst the methods discussed above, the moment-matching methods of Satterthwaite-Welch, Hall-Buckley-Eagleson, Wood and Lindsay-Pilla-Basak are the only options that meet these criteria and are described in Sect. 3 below. The first three of these methods only require a single evaluation of a particular cdf and the storage of a fixed number of parameters that can be easily sequentially updated. The Lindsay-Pilla-Basak method is more computationally intensive, but has the potential to give more accurate results by matching higher-order moments. There are other approximate methods (e.g. Solomon and Stephens 1977) besides these four, but they all have shortcomings (e.g. require too much memory, too expensive to compute) that would render them unsuitable for streaming data applications.
Our motivating application for computing F Q N is as part of a sequential change detector for the variance of a process; see Bodenham (2014, Chap. 8) for methodological background, and Ye et al. (2002) for an application in computer network security. Suppose we are interested in making inference on the sequence z 1 , z 2 , . . . , z N which are observations generated from random variables Z 1 , Z 2 , . . . , Z N , and the weighted variance is defined as where c = (c 1 , c 2 , . . . , c N ) are some weights, andZ is the (possibly weighted) mean of Z 1 , . . . , Z N . If the Z i are i.i.d normal, then it can be shown that V c,N is distributed as some Q N . This formulation is similar to the exponentially weighted moving variance described in MacGregor and Harris (1993). In a streaming data scenario, it would be infeasible to use a method such as Imhof's which requires the storage of the whole vector c, particularly when N becomes large. In sequential change detection, N increases until a change is detected. The size of the change would then depend on the application; in cybersecurity problems of interest to us, we expect N to be between 100 and 1000. Streaming data algorithms need to have low and fixed memory requirements and be computationally inexpensive.

Efficient approximate moment-matching methods
As the name suggests, these methods involve matching the moments of Q N to those of another distribution, and using that distribution's cdf to approximate F Q N . In order to do this, the moments of Q N need to be computed. However, instead of computing the moments directly, it is easier to first compute the cumulants of Q N and then obtain the moments from the cumulants. In fact, the first three methods described below directly use the computed cumulants, and do not require computation of the moments.

Computing cumulants and moments
The cumulants of Q N , a weighted sum of i.i.d. χ 2 1 random variables as in Eq. (1), are denoted by κ r (Q N ) and can be computed using the formula where d = (d 1 , d 2 , . . . , d N ) are the weighting coefficients. This can easily be shown using the properties of cumulants and recalling that for a χ 2 1 random variable X , κ r (X ) = 2 r −1 (r − 1)! [e.g. Box (1954)]. In a sequential context, when Q N becomes Q N +1 , the cumulants can be easily updated by For the remainder of this chapter, we shall only be concerned with Q N , and so shall write κ r = κ r (Q N ). The moments of Q N , denoted m r = m r (Q N ), can be computed from the cumulants using m 1 = κ 1 and Since the first three methods described below only require the first two or three cumulants of Q N , these are explicitly provided here:

Satterthwaite-Welch approximation
Equating the first two moments of Q N with a Γ ( k, θ) variable yields If we use F Γ (k,θ) to denote the cdf of a Γ (k, θ) distribution, then the Satterthwaite-Welch approximation uses F Γ ( k, θ) to approximate F Q N . In the references [e.g. Box (1954)], the Γ (k, θ) distribution is often written as a scaled χ 2 1 distribution.

Hall-Buckley-Eagleson approximation
We provide a brief outline of the method which is fully described in Buckley and Eagleson (1988). First, Q N is used to denote Q N normalised as in Second, if ν is defined as and

Wood F approximation
Wood's F method (Wood 1989) matches the first three moments of Q N with another distribution that has a probability density function of the form where is the beta function. Although in Wood (1989) it is referred to as an F distribution, the density in Eq. (12) can be better described as that of a G3F or corrected F distribution (Pham-Gia and Duong 1989; Johnson et al. 1995). The parameters α 1 , α 2 , β can be defined in terms of the cumulants κ 1 , κ 2 , κ 3 computed in Eq. (4) above (e.g. using Gröbner bases): It is noted in Wood (1989) that if X is distributed according to the density in Eq. (12), then where F(2α 1 , 2α 2 ) is a standard F-distribution with parameters 2α 1 and 2α 2 . Therefore, if Y ∼ Q N , and y is an This approximation can be used as long as both r 1 , r 2 > 0, which is guaranteed in many cases (Wood 1989). When either r 1 = 0 or r 2 = 0 (it is proved in Wood (1989) that neither can be negative), then Wood (1989) recommends using either the Satterthwaite-Welch approximation, or another two-moment approximation.

Lindsay-Pilla-Basak approximation
The method described in Lindsay et al. (2000) approximates where each π i ≥ 0 and i π i = 1, and the 2n +1 parameters k, θ 1 , θ 2 , . . . , θ n , π 1 , π 2 , . . . , π n are to be determined. These parameters are computed by following a sequence of steps that make use of results concerning moment matrices (Uspensky 1937, Appendix II). The sequence in Lindsay et al. (2000) is complicated, so we extract the main steps here (without proofs). The first step is to compute the first 2n cumulants κ 1 , κ 2 , . . . , κ 2n of Q N using Eq. (4), and then use the recursive formula in Eq. (6) to compute the first 2n moments m 1 , m 2 , . . . , m 2n of Q N . The second step is to define, for a variable α, the functions δ r (α) as and δ 0 (α) = 1. These functions are then used to create the (r + 1) × (r + 1) pseudo-moment matrices r (α), defined as For example, The third step is to find certain roots λ 1 , λ 2 , . . . λ n such that For r = 1, there is a unique positive root λ 1 = m 2 /(m 2 1 ) − 1. For r > 1, one can use a bisection method (e.g. Everitt 2012) to solve for the root λ r ∈ [0, λ r −1 ) of the equation det r (α) = 0. Eventually, λ n is obtained. The fourth step is to define the matrix M n ( λ n , t), Note that M n ( λ n , t) is the same as n ( λ n ) but with the last column replaced by (1, t, . . . , t n ) . This matrix is used to compute the nth degree polynomial S n (λ, t), where for some c j ∈ R and j = 0, 1, . . . , n. In order to obtain the value of the coefficient c j , one can replace the last column of M n ( λ n , t) (the powers of t), with the basis vector e j+1 (the ( j + 1)th component equals one, all others are zero), and compute the determinant of this modified matrix. With the coefficients computed, the n roots of S n (λ, t) = 0, denoted μ 1 , μ 2 , . . . , μ n , can be found [the roots are real and distinct (Uspensky 1937, Appendix II.4)]. The fifth step is to use these roots μ i to solve the system of linear equations . . .
to compute the mixture proportions π 1 , π 2 , . . . , π n . Since the matrix on the left of Eq. (25) is a Vandermonde matrix, it is non-singular (Macon and Spitzbart 1958), and so this system of linear equations has a unique solution. Finally, we define k = ( λ n ) −1 and θ i = λ n · μ i , for i = 1, 2, . . . , n, and now can compute the approximate cdf F Q N in Eq. (17). Note that the Lindsay-Pilla-Basak method agrees with the Satterthwaite-Welch method for n = 1. It should be remarked that Robbins and Pitman (1949) also attempt to obtain an approximation using a method of mixtures, but by computing the characteristic function rather than using the method of moments.

Sequential implementation
As described in Sect. 2, one might wish to compute F Q N (x) and then soon afterwards compute Note that x and x may be different values. This can be done easily and efficiently using one of the four moment-matching methods described above. When comput- where the value of depends on the method we are using (e.g. for Hall-Buckley-Eagleson, = 3). Now, one can simply use the new coefficient d N +1 and Eq. (5) to update κ r (Q N ) to κ r (Q N +1 ), for r = 1, 2, . . . . These updated cumulants, together with x , are all that is needed to compute F Q N +1 (x ). Note that this method only requires the storage of the cumulants, regardless of the value of N , which makes this method suitable for a streaming data context.

Evaluation of approximate methods for computing F Q N in the literature
In previous work on approximations for computing the cdf F Q N of weighted sums of chi-squared random variables Q N (Imhof 1961;Solomon and Stephens 1977;Wood 1989;Lindsay et al. 2000;Castaño-Martínez and López-Blázquez 2005), it was common to estimate the performance of an approximate method by demonstrating its accuracy for a selected sample of M distributions and d k = (d 1,k , d 2,k , . . . , d N ,k ). Recall that the cdf of a random variable X is defined by In this article, values x in the domain of the cdf F X will be called quantile values, and values F X (x) will be called probability values. For each Q N ,d k , the quantile values x j,k are found such that, for k = 1, 2, . . . , M, for a specific set of probability values p j . Then a table of errors j,k , where, for k = 1, 2, . . . , M, is presented for one or more approximate methods, where G is the cdf produced by the approximate method. According to the literature, the method with the smallest set of errors is then considered to be the best approximate method. This may seem to be a reasonable approach, but the execution in previous works leaves something to be desired. In Imhof (1961), Solomon and Stephens (1977), Wood (1989), Lindsay et al. (2000), Castaño-Martínez and López-Blázquez (2005), each analysis only considers a selection of between M = 8 and M = 18 distributions Q N for a selected set of coefficients and number of terms. Results established for an approximation procedure based on the analysis of such a small selection should be viewed with caution. So, while previous works may have established the accuracy for the particular selections considered, those results cannot reasonably be assumed to hold for all possible Q N . Moreover, previous works only considered Q N with fewer than N = 10 terms, so it is natural to wonder how approximate methods perform for distributions Q N with significantly larger N . This is particularly relevant in the context of streaming data problems.
There is a possible explanation for why previous works only consider a limited selection of distributions Q N in their analyses. When these approximate methods were first considered in the 1950s and 1960s (e.g. Box 1954;Imhof 1961), calculating the probability values p j may have been difficult, especially with computing in its infancy. Therefore, only a limited table of results was produced. When later methods in the 1970s and 1980s (e.g. Solomon and Stephens 1977;Wood 1989) were developed, it would have been natural to use the performance analysis of earlier methods as the benchmark, and so a table of errors j,k was again compiled for a small (in some cases the same) sample of distributions. Unfortunately, this method of evaluating performance has continued unchanged (e.g. Lindsay et al. 2000;Castaño-Martínez and López-Blázquez 2005), even though computers that are able to complete a much more thorough analysis are now readily available. In Sect. 5, we outline such an analysis, which will seem natural following the discussion in this section.
It should be mentioned that while we shall use Farebrother's method in combination with a bisection procedure (e.g. Everitt 2012) to compute the exact quantile values [i.e. Eq. (28)] in Sect. 6, it was not indicated in previous works how the exact quantile values were obtained for performance calculations.

A new method for evaluating the performance of an approximate method for a cdf of a weighted sum of random variables
This section discusses the issue of evaluating the performance of approximation methods for the cdf of a weighted sum of random variables. This procedure is then used in Sect. 6 to analyse the performance of approximate methods for the cdf of a weighted sum of chi-squared random variables. In this section, R N is a weighted sum of i.i.d. unspecified random variables (not necessarily chi-squared as Q N ). It is assumed that a method exists for computing the true probability value F R N (x) for quantile value x, to arbitrary accuracy. However, the method may be too computationally or memory intensive for routine application.

Performance of an approximate method for a particular distribution R N,d
Suppose a method provides approximate probability values G(x) for a weighted sum of random variables R N . Suppose further that we wish to determine how close G is to the true cdf F R N , for a particular distribution R N ,d with weights d =  (d 1 , d 2 , . . . , d N ). For a set of probability values suppose that the "exact" quantile values can be computed to an arbitrary precision, perhaps at a practically unacceptable computational cost, so that In this case, we shall say that the quantiles are accurate to precision ξ , when we mean that the true cdf will evaluate the quantile to within ξ of the corresponding probability value. The errors of the approximate method G, denoted by j , are then defined as The smaller the j , the better that G approximates F R N for the probability values p j . By a simple application of the triangle inequality, is obtained. Therefore, if the x j can be computed to ensure ξ j for all j, it is then only necessary to look at the values |G(x j ) − p j | to obtain a good approximation for j .

Estimating the accuracy of an approximate method for R N , for a particular N
The first step to more comprehensively evaluating the performance of an approximate method for distributions with N terms is to randomly generate a large sample of M coefficient vectors d k = (d 1,k , d 2,k , . . . , d N ,k ), where for some distribution D, so that F R N ,d k is the cdf of for some distribution Y . The next step is to select a wide range of probability values { p 1 , p 2 , . . . , p L }, and then to compute the quantile values so that for some precision ξ , with ξ 1, Finally, the errors j,k are computed as for j = 1, 2, . . . , L and k = 1, 2, . . . , M. The set of errors for probability value p j is defined as E j = j,k |k = 1, 2, . . . , M , j = 1, 2, . . . , L . (40) While it would now be easy to compute max E j and declare this to be a reasonable upper bound for the error when computing p j , provided that M is large, the following procedure is preferable because it establishes a probabilistic result. Define¯ j to be the sample mean, s 2 j the sample variance, and q 2 j the scaled sample variance of E j by the equations: Suppose that * j is the error for F R N ,d * , with coefficient vector d * generated as in Eq. (35). If we assume that the error values in E j are i.i.d. according to some distribution, then Chebyshev's inequality with the sample mean and variance Saw et al. (1984) gives us, for any δ > 0, If we set the the right-hand side of Eq. (44) to be then Eq. (44) implies Then¯ j + δq j provides an upper bound for 100(1 − α δ,M )% of all possible errors obtained when computing p j using the approximate method. In other words, the probability that the error exceeds the upper bound is less than α δ,M . For example, when δ = 10 and M = 10,000 then α δ,M ≈ 0.01, or δ = 32 and M = 10,000 gives α δ,M ≈ 0.001, and so then¯ j + δq j provides an upper bound for 99.9% for all errors. The same procedure could be followed to obtain a bound for the error of computing p j for every p j ∈ { p 1 , p 2 , . . . , p L }, and so an estimate of the error for an approximate method of computing probability values for distributions Q N is obtained, for a particular N .
The assumption that the errors in E j are i.i.d. may seem restrictive, but in fact the errors need only be weakly exchangeable. Finally, although Saw et al. (1984) give a slightly sharper bound for the inequality in Eq. (44), its expression is far more complicated and does not significantly change the bound for our purposes here.

Results
For the purposes of discussion below, let us define the lower tail to be the probability values in P L and the upper tail to be probability values in P U . Values in P M will be referred to as middle probability values. Assuming that the coefficients are sampled from U (0, 1) is not particularly restrictive; if a particular application uses coefficients that are known to be bounded, they can be rescaled to the range (0, 1). Farebrother's method is used to ensure that the quantiles are accurate to ξ = 10 −8 as in Eq. (38). Imhof's method could also have been used, but the implementation of Farebrother's method in the R package CompQuadForm (Lafaye de Micheaux 2011) appears to allow a greater precision to be specified. The analysis is then performed using δ = 32 to obtain an upper bound with confidence α δ,M ≈ 0.001 (see Eq. (45)). The accuracy of each of the four momentmatching methods in Sect. 3 is computed for all p j , and the methods are compared side by side in Sect. 6.1. The Lindsay-Pilla-Basak method is computed for n = 4 (that is for the first eight moments), and so will be abbreviated to LPB4. In Sect. 6.4, we then investigate the relative speeds of each method. Note that none of the sampled coefficient vectors d k yielded degenerate cases (as mentioned in Sect. 3.4) for the Wood F approximation.

Accuracy
The  Figure 1 illustrates several points. The first feature of interest is that the methods generally increase in accuracy as N increases. There are a couple of exceptions (e.g. p j = 0.999 for the LPB4), but any decrease is minor. This seems to suggest a trend which would continue for N ≥ 100 (indeed, similar figures showing results for N = 200, 500 and 1000 confirm this). Following this observation, if method A has number of digits of accuracy a for number of terms N , we shall say that method A is accurate to a decimal places for N ≥ N . As far as we are aware, this observation that the accuracy of these approximate methods generally increases, as the value of N increases, has not been noted before and is not apparent or implied from the construction of the methods. As already mentioned, previous analyses only focused on distributions Q N for a limited range of N .
If the results shown in Fig. 1 for each individual method are now examined, it can be seen that SW is accurate in the upper and lower tails to at least two decimal places for N ≥ 100.  The HBE method is accurate to two decimal places for all p j for N ≥ 50 and to three places for almost all values in the upper and lower tails for N ≥ 100. The WF method is also accurate to two decimal places for p j for N ≥ 50 and is accurate in the upper tail to three digits for N ≥ 50. The LPB4 method is accurate to four decimal places for almost all prob- Fig. 3 Error of the Normal approximation, compared to the Satterthwaite-Welch approximation, grouped by method ability values (only exceptions are a few middle probability values) for N ≥ 50 and has close to five digits of accuracy for the upper and lower tails for N ≥ 100. Note that Fig. 1 is meant to illustrate the general behaviour for each method across a range of probability values, as N increases. In the supplementary material, this figure has been split into three separate figures, displaying the upper, middle and lower probability values, for readers who may be interested in the behaviour for a particular probability value. Figure 2 shows that over the different probability values, SW is the least accurate, while LPB4 is clearly the most accurate, and WF and HBE appear to be essentially matched, although for most probability values WF has a slightly better accuracy than HBE (one exception is for p j = 0.975 and N = 50).
Note that if Imhof and Farebrother's methods were included in Figs. 1 and 2, since they are essentially exact (they will iterate until the desired accuracy is achieved), the result would be horizontal lines at the level of the accuracy specified.
One reviewer raised the question of how these methods perform for very small probability values. An investigation into the accuracy of these methods for small probability values in the set {10 −4 , 10 −5 , . . . , 10 −10 } is included in the supplementary material, which shows that the Wood F and Lindsay-Pilla-Basak methods perform well for probability values in this range, but that the Hall-Buckley-Eagleson method should not be used in this case.
Another reviewer raised the question of how these methods perform for coefficients that are not U (0, 1)-distributed or are not i.i.d.. A section in the supplementary material shows similar performance to that in Fig. 1 for coefficients that are Beta(2, 5)-distributed, for coefficients that are sampled from a mixture of distributions, and for coefficients that are highly correlated. These results indicate that the actual distribution of the coefficients is not too important when considering the results in Fig. 1. Finally, other sections in the supplementary material show similar results for which the variables are χ 2 (ν) for ν > 1, rather than χ 2 (1), and that the accuracy of the methods increases when N = 200, 500, 1000.

Comparison to the normal approximation
Although the normal approximation is not considered to be as good as the four approximations considered above, it is interesting to investigate how it compares to SW, the simplest of the approximations above.
The normal approximation is computed in a similar manner to SW. Equating the first two moments of Q N with a N( μ, σ 2 ) variable yields following the definition of the cumulants κ 1 and κ 2 in Sect. 3.1. Then F N( μ, σ 2 ) is used to approximate F Q N . Figure 3 shows that SW appears to be one decimal place more accurate than the normal approximation. The only exception is for p j = 0.999, where the two methods appear to have similar accuracy. Even though both methods are two-moment approximations, and the computational complexity is virtually the same, SW's use of a Gamma cdf provides a significant increase in accuracy over the normal approximation.

Accuracy for small number of terms N
It is worth investigating the accuracy of these methods for the cases where N ∈ {2, 3, ..., 10}. The results of this investigation are shown in Fig. 4 and show that SW, HBE and WF will generally give between 0 and 2 digits of accuracy, while LPB4 generally gives at least 2 digits of accuracy. These results suggest that when N < 10, these methods should be used with caution. Note that for N = 2, 3 there are choices of coefficient vector d k which result in the LPB4 method not being able to provide an approximation [fails to find roots λ r for Eq. (22)], so values for N = 2, 3 for LPB4 are omitted. If one were only interested in computing the cdf for a fixed, small N , then, as one of the reviewers has suggested, Imhof's method should be used. However, if the number of terms N were increasing, as in a change detection scenario (see Sect. 2), it would be better to use one of the momentmatching methods for all N . Table 1 shows that while the SW, HBE and WF methods have similar speeds (of the same order), LPB4 is significantly slower. This could be due to the iterative methods needed in steps 3 and 4 of the algorithm (as described in Sect. 3.5) and the matrix algebra in several steps. Besides the matrix operations, the LPB method needs to employ root-finding algorithms (which can be very efficient, but are still iterative). For comparison purposes, the speeds of the normal approximation, Imhof's method and Farebrother's method have also been included. The normal approximation is slightly faster than SW, but is much less accurate. Surprisingly, Imhof's method is faster than LPB4, but is still over 40 times slower than HBE. LPB4 is over 300 times slower than HBE. Farebrother's method is significantly slower than any of the other methods, but it is unclear if this is due to a few problematic cases, or if this is a general property of the algorithm. However, the table shows its performance over 10,000 samples, which gives an indication of its average behaviour. The four algorithms (SW, HBE, WF and LPB) and the normal approximation were written in R, while Imhof's method and Farebrother's are implemented in C++ in the R package CompQuadForm (Lafaye de Micheaux 2011). Note that the implementation in C++, a compiled language, may explain why Imhof's method is faster than LPB4. The speed test was done on an Apple iMac with an Intel Core i5 (3.2 GHz) processor (4 cores) and 8 GB of RAM.

Conclusion
While Imhof's method is essentially exact, it is not suitable for a streaming data scenario, where it is necessary for algorithms to (a) not store all the coefficients of Q N , and (b) have efficient computation. In such situations, moment-matching methods such as the four described in Sect. 3 may be very useful.
Choosing between these methods is not a simple matter of choosing the most accurate. One also needs to consider the speed of computation, and, to a lesser extent, the ease of implementation. While Figs. 1 and 2 show the Lindsay-Pilla-Basak method to be extremely accurate, it is also significantly slower to compute (see Table 1) and laborious to implement (Sect. 3.5). If it is not necessary to have four decimal places of accuracy, other methods could be used.
Of the remaining three methods, the Hall-Buckley-Eagleson method is perhaps the best alternative. It is one decimal place more accurate in the tails than the Satterthwaite-Welch method, yet is only marginally slower (see Sect. 6.1 and Table 1), and is essentially as accurate as the Wood F method, without needing to worry about degenerate cases (see Sects. 6.1 and 3.4). For this reason, the Hall-Buckley-Eagleson method is recommended for most practitioners.
This recommendation is based on the observation, revealed by Figs. 1 and 2 and not previously described in the literature, that the accuracy of the four moment-matching methods generally increases as the number of terms N increases.
However, as described in Sect. 6.1 and shown in the supplementary material, for very small probability values, either the Wood F or the Lindsay-Pilla-Basak method should be used.
Furthermore, Sect. 5 provides a new statistical framework for evaluating the accuracy of an approximate method for computing F R N , the cdf of a weighted sum of random variables R N (for any distribution).