 # Bootstrapping Through Discrete Convolutional Methods

Bootstrapping was designed to randomly resample data from a fixed sample using Monte Carlo techniques. However, the original sample itself defines a discrete distribution. Convolutional methods are well suited for discrete distributions, and we show the advantages of utilizing these techniques for bootstrapping. The discrete convolutional approach can provide exact numerical solutions for bootstrap quantities, or at least mathematical error bounds. In contrast, Monte Carlo bootstrap methods can only provide confidence intervals which converge slowly. Additionally, for some problems the computation time of the convolutional approach can be dramatically less than that of Monte Carlo resampling. This article provides several examples of bootstrapping using the proposed convolutional technique and compares the results to those of the Monte Carlo bootstrap, and to those of the competing saddlepoint method.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Bootstrapping is a resampling technique that relies on taking random samples with replacement from a data set. Bootstrapping techniques (Efron, 1979; Efron and Tibshirani, 1994) provide researchers with an increased capacity to draw statistical inference from a single sample, whether or not parametric assumptions are made. The theoretical bootstrap distribution of a statistic exists, but as sample size (n) increases, the number of possible bootstrap samples grows quickly with possibilities. Monte Carlo methods are standard practice for approximating the bootstrap distribution of a statistic.

While the Monte Carlo bootstrap is extremely versatile, as with all procedures, it has some limitations. Monte Carlo bootstrap methods are relatively straightforward to implement but can be computationally intense, and produce stochastic error bounds which narrow at the relatively slow rate of . This article introduces an alternative method for bootstrapping quantities using convolutional methods which, when short computation times and high levels of accuracy are valued, often proves to be better suited than the Monte Carlo approach and the competing saddlepoint method. The convolutional method is computationally fast and provides either exact values or more precise approximations.

For a simple illustration, consider the sample mean, , from a sample {} for . The sample implicitly defines a discrete distribution with the support

with each atom being equally probable (assuming no ties). If for instance, one is interested in finding the bootstrap distribution of the mean, it is primarily characterised by self-convolving a modified version of that distribution

times. Since the distribution of the bootstrap mean is defined by the sum of discrete random variables with compact support, it also has a discrete distribution with compact support. This simple scenario outlines the basics of the method described in this article.

The numerical techniques for computing convolutions are numerous, however in this work we focus on using the discrete Fourier transform (DFT). The formulae for computing both the DFT and its inverse have accompanying algorithms in the fast Fourier transform (FFT). These algorithms are readily available on most, if not all, computational platforms. The FFT scales extremely well according to complexity (Cooley and Tukey, 1965). In some instances, the method we propose can compute the bootstrap distribution exactly. However, in the more general case, when discretization error is introduced, mathematical bounds can be obtained for quantities of interest.

We provide an outline of the remainder of the paper. First, a review of the relevant literature is presented. Once the literature has been discussed, the methodology is explained. The next section focuses on a comparison between the discrete convolutional and Monte Carlo methods. The comparison is followed by a section containing five data applications of the discrete convolution method. The examples also serve as an additional comparison between computational approaches. Finally, in the conclusion we summarize the main points of this paper.

### 1.1 Literature Review

Exploring alternatives for Monte Carlo bootstrapping techniques is not uncommon. There have been a few lines of research that have produced excellent results. In this section we list some of these works and their relevant findings.

The most prominent alternative approach to approximate bootstrap distributions is the use of the saddlepoint method. In those endeavors the saddlepoint approximation inverts the moment generating function of a random variable to obtain an accurate estimate for the probability density function. For example,

Davison and Hinkley (1988) developed methods to estimate the bootstrap distribution of a statistic using saddlepoint approximations. Other work along these lines includes Butler and Bronson (2002)

, which bootstrapped survival curves using a saddlepoint approximation. This same paper focused on using the saddlepoint method to estimate the mean and variance of first passage times (FPT). In an example studying the progression of dementia, the saddlepoint method was generally able to predict values closer to the truth than the Monte Carlo bootstrap.

More recently, the saddlepoint method was used to approximate first passage quantiles of a stochastic process. For example,

Balakrishnan and Qin (2019) developed a nonparametric method to approximate FPT quantiles of a degradation process. This methodology was further improved in Palayangoda et al. (2020) which increased the accuracy and proposed a method to handle unequally spaced data. Given the popularity of the saddlepoint method, it will be used as one of the main comparisons for the examples in this article and we will compare our proposed methods directly with results from both of these works as well as others.

The method proposed in this article is based on the discrete Fourier transform and its celebrated computational implementation, the fast Fourier transform. The convolutions inherent in a bootstrap distribution are computed in the Fourier transform space. Previous work has applied the FFT to provide exact p-values (Beyene, 2001). FFT methods have been used to calculate the distribution of a random variable, as well as find exact hypothesis test p-values for general linear models.

Additional research has calculated the distribution of convolutions with the FFT. Kern et al. (2003) used the FFT to approximate a two-dimensional distribution that characterizes animal locations. Although the estimated distributions are different from those approached in this paper, use of the FFT was shown to drastically decrease computation time. Additionally, Warr (2014); Warr and Wight (2020) provided mathematical bounds on the convolutions of random variables using the FFT. We use these established methods to accurately calculate quantities from bootstrap distributions.

Recent research has also developed the algorithmic construction of bootstrap confidence intervals (Efron and Narasimhan, 2020). Although the goal of that research is not to revolutionize the estimation method of a bootstrap distribution, it demonstrates that bootstrapping methods continue to improve and are becoming an ever more viable option for inference. Furthermore, Efron and Narasimhan state that the bootstrap confidence interval often enjoys better accuracy over intervals constructed using large-sample normal approximations.

## 2 Methodology

The typical Monte Carlo bootstrap is performed by randomly resampling from the data. In this section we propose an alternative method that computes the bootstrap distribution directly.

### 2.1 Bootstrap Distributions

To establish notation we define a distribution, , such that we have mutually independent random variables for . Once the random variables are drawn from their distribution, they are no longer random, and we denote the observed sample as or just . Now a discrete probability measure can be defined from the data such that

 Y∼G , where G=1nn∑i=1δxi.

We define to be the Dirac measure at . The sample now defines a discrete distribution.

The random variable itself is typically not of interest, however, we are often interested in statistics which are functions of . Thus we denote a generic bootstrap statistic of as . Our proposed method is intended for cases where the bootstrapped statistic can be written as the sum of independent random variables.

### 2.2 Convolutions using the Discrete Fourier Transform

The sum of independent random variables is an example of a convolution. Similar to using moment generating functions, the discrete Fourier transform provides a readily available technique to calculate the convolution defined by the sum of independent random variables.

The discrete Fourier transform is the primary tool in our method to calculate the distributions of bootstrap statistics. To use the DFT, an equally spaced grid of length is defined on the support of the random variable . This grid must be defined such that it adequately covers the support of and . We denote the grid on the support to be . Then the DFT for is the sequence defined by:

 ~fY,k=N−1∑j=0P(Y=sj)e−i2πNkj

for .

If were to be convolved with a second random variable, , the DFT for is the sequence defined by , so long as and are independent.

In order for this convolutional method to be exact, must precisely align with the supports of , and . Here we note that the grid and the actual supports of those random variables may not align. If this is the case, an approximate distribution for can be found and the error bounds of the approximation can be obtained. It is also important to point out that random variables defined by a sample, such as , have discrete and compact support. The discreteness property aides in finding an appropriate distance between points in , and the compactness property allows for straightforward selection of a starting and end point for .

The basic approach for bounding the error of an approximated convolution is: if all exact support points (not contained on the grid ) of the random variables are rounded down to the next lower grid point, and the DFT convolution method is used, the resulting CDF will be an upper bound for the distribution. Likewise, if all exact support points of the random variables are rounded up to the next greater grid point, the result is a lower bound for the CDF. Using this approach, it is possible to mathematically bound the true convolved distribution in cases where it isn’t practical to ensure the support of the random variables and the grid completely agree.

One desirable feature of the DFT is the ease in calculating both the forward and inverse DFT using a fast Fourier transform algorithm. Using the FFT and its inverse, it becomes a simple calculation to find the DFT, convolve in the Fourier domain, and then use the inverse transform to find the CDF of the bootstrap statistic.

The process for finding the bootstrap distribution of a statistic (which is defined by a sum of independent random variables) can be described in a few simple steps. First, the random variables to be convolved should be identified. These will be some function of the data. For these statistics, the convolved random variables should be independent, but don’t necessarily need to be identically distributed. Next, a grid, , should be defined. This grid needs to be large enough to contain the full support of and the bootstrap statistic, . The grid consists of equally spaced points. Note that the number of grid points will determine the accuracy of this method. In general, more grid points will produce more accuracy. When the support of the distribution and the grid agree, the resulting distribution will be the exact distribution of the bootstrap statistic. Also note that some FFT algorithms are more efficient if the number of the grid points in the support is a power of two. Third, the probability mass functions for the random variables to be convolved must be defined. The probability mass functions will need to agree with or be fit to the the support grid, . Next, the DFT for each random variable is found (utilizing the FFT). The DFTs can then be multiplied together to perform the convolution. Finally, by inverting the resultant DFT, the pmf of the bootstrap distribution is obtained. This process will be demonstrated several times in the examples and applications.

### 2.3 The Bootstrap Distribution of the Mean

The sample mean is one of the most commonly bootstrapped statistics. The distribution of the bootstrap mean can be defined as follows:

 ¯Y=1nn∑i=1Yi, for Yi %iid∼ G. (1)

The distribution of provides an estimate for the population mean (with reasonable stochastic uncertainty) when the distribution has finite variance (Knight, 1989).

One can certainly obtain Monte Carlo samples from the distribution of , but the discreteness of invites the use of discrete methods to obtain the distribution of , which reflects our uncertainty regarding the true population mean. This can be done using the convolution defined in Equation 1 and the discrete Fourier transform. A consequence of the discrete nature of is that the distribution of and various other random quantities also have discrete distributions.

For this example let . The shifting of ensures that the support is simple to work with. In order to employ the DFT to find the bootstrap mean’s distribution, we define an equally spaced grid, , of length . When choosing , we attempt to include all the values and define it on the interval . The ’s define a probability measure, , on such that each contributes mass of . Many grid points, , will have a mass of 0. Let the random variable and denote the DFT sequence of for . The DFT sequence of , where the ’s are independent draws from , is:

 ~fV,k=(~fZ1,k)n. (2)

The pmf for the random variable is simply found by invoking the inverse discrete Fourier transform on the sequence defined in Equation 2. A common definition of the inverse DFT is:

 P(V=sj)=1NN−1∑k=0~fV,ke−i2πNkj (3)

for . Note, that after employing the inverse FFT, some implementations omit the term. This is certainly the case in R (R Core Team, 2021), and the result must be divided by afterwards. Once shifted by the pmf of becomes the bootstrap distribution of the mean, or in our more standard notation the distribution of .

The primary issue with employing this convolutional method is that it requires an evenly spaced grid on the support which includes all values . This is often not tenable in practice. however, it is possible to create a grid that nearly contains each value and the mass assigned to that value can be reassigned to the nearest grid point. While some precision will be lost, this process provides a readily obtainable approximation. We also show how to bound the error of that approximation.

### 2.4 Illustrative Examples

This first example investigates an extremely simple case. However, the example should provide some intuition into how the convolutional method works.

Consider the sample . In this instance, due to the small sample size, it is easy to find the exact bootstrap distribution of by enumerating the sample space. The exact distribution is displayed in Figure 1. Figure 1: In the case of just 4 observations, it is relatively simple to find the exact bootstrap distribution of the sample mean.

To apply the discrete convolutional method described in this paper, the support needs to be broken into a grid of equidistant support points. For this example, the grid was defined as the sequence of numbers from 0 to 8 increasing in increments of 0.25. This support was chosen as it includes each and covers the entire support of . The results are exact (up to numerical precision).

After dividing the observations by the sample size of 4 there are now four observed contributions that match with four of the predetermined grid points (0.25, 1.00, 1.50 and 2.00). A probability sequence of the same length as the support should now be generated. Where an observation matches a support point, the value in the new sequence will be 1/4, while all other values will be 0. That is, the term 1/4 should be found in the second, fifth, seventh and ninth positions. This characterizes the pmf of the discrete random variable .

Now properties of the discrete Fourier transform can be used to find the bootstrap distribution of the sample mean. This sequence of probabilities should be transformed using the DFT. Next the sequence will be raised to the power, in this case the fourth power, to account for the convolution. This results in the DFT for the bootstrapped sample mean. Finally, the sequence can be un-transformed, using the inverse FFT, and it provides the exact pmf of .

This is one of the simplest applications of the convolutional method. An additional step in complexity is added when there is no obvious way to create an equally-spaced grid that contains the entire support of the statistic of interest. Now consider changing the previous example by replacing an observation such that .

An intuitive approach would be to create a grid and then round observations to the closest grid points. However, we advocate for a slightly different approach, which allows the error of the approximation to be mathematically bounded. The process for finding an approximation is similar to the exact process shown above, but the process is repeated twice, once for the upper bound and once for the lower bound.

To create a lower bound on the cumulative distribution function (CDF) of interest, rather than rounding the data, each data point can be shifted up to the nearest higher grid point. Likewise, finding the upper bound shifts the mass of each observation to the nearest lower grid point.

The theory for the process of bounding a distribution through convolutional methods was presented in Warr and Wight (2020). By computing the difference of the upper and lower bounds of the CDF at any of the grid points on the support, a mathematical bound is placed on the distribution of .

The width of these bounding intervals can be changed by increasing or decreasing the number of grid points. With more support points the interval widths will be narrower.

Referring back to the example; to create the bounds, the grid was taken as a sequence from 0 to 9 with a difference of 0.1 between adjacent terms. Once the grid was defined, the data were fit to the grid as described above and the convolutional method was employed to place bounds on the distribution. Figure 2: The CDF for the bootstrap distribution of the mean can be bounded. The upper bound is given by the solid line, while the dashed line is the lower bound. The true CDF is contained between the upper and lower bounds.

As shown in Figure 2, even when it isn’t feasible to find the exact bootstrap distribution, the distribution can be bounded. It would be simple enough to increase the number of grid points which would result in a more accurate approximation of ’s distribution.

In addition to the sample mean, the convolutional method is effective in bootstrapping other quantities. In this section we discuss an approach to find quantiles of the bootstrap distribution of a product’s failure time using degradation data. Some of the applications presented later in this article demonstrate this in detail.

In the simplest case, consider a collection of items that experience a random amount of degradation , for , which occurs during a fixed time period . That is, if the system is inspected in time increments of , the increase in cumulative degradation since the last inspection is drawn from the distribution, . In this case we consider the to be independent and identically distributed random variables. Once observed, the define a discrete random measure and we assign the random variable to that distribution. If the total amount of degradation meets or exceeds some predefined threshold, , then we consider an item to have failed.

Thus when using Monte Carlo sampling to bootstrap degradation data for failure times, we sample independent until their sum meets or exceeds . Suppose we sample ’s, such that and . Then the sampled item’s failure time, , is calculated to be

 Z=d∗[(K−1)+(T−K−1∑i=1Yi)/YK].

Notice that is a random variable. We are defining the distribution as the distribution of bootstrapped failure times. Since the degradation is only observed at time increments of , the definition of

uses a linear interpolation to generate a failure time.

A quantile can be estimated directly from one bootstrap sample, or multiple samples can be obtained which provides some Monte Carlo uncertainty of the estimate. Obtaining multiple samples can be computationally intense. Suppose we estimate the distribution of failure times with Monte Carlo samples. The resulting draws can be used to estimate any quantile, however in order to assess uncertainty, the process must be repeated many times to produce multiple draws for the quantile of interest. One issue that makes the Monte Carlo approach less appealing when estimating quantiles is that the Monte Carlo error increases quite dramatically when finding quantiles far from 0.5. Our proposed procedure to estimate quantiles is much faster and more accurate than the Monte Carlo approach.

To calculate this bootstrapped failure time distribution using the convolutional method we find the CDF directly. For a given time we find a value for . To do this, first find the integer such that and . Next, define a grid, , of length , where and . Again, controls the accuracy of this approach, where larger is preferable. Next, obtain the DFT of and , where . Now the approximated DFT sequence of is

 ~fZ,k=~faY,k(~fY,k)k−1.

The final steps are to invert the DFT sequence of and then sum the probability mass on the support grid which is less than . Quantiles can found by trying different values of .

## 3 Comparisons

One advantage of using the convolutional method is that the computation time can be much shorter when compared to traditional Monte Carlo methods. The time needed to generate bootstrap intervals using the Monte Carlo method is based on the number of bootstrap samples. In contrast, the computation time for the discrete convolutional method is based on the length of . When a small number of bootstrap samples are taken using the Monte Carlo approach, there is little difference in computation time for the two methods. However, in some cases a large number of Monte Carlo samples are required and here the convolutional method will be faster with more accuracy.

The accuracy of the convolutional method is determined by the agreement of and the support of the sample. In cases where the support of the sample is given to a fixed precision (e.g., rounded), the user defined grid, , can often be defined to include all possible bootstrap outcomes. In these instances, increasing the number of support points on the grid will not improve the results (since they are exact). When a user defined grid cannot include all possible bootstrap samples, mathematical bounds can be placed on the bootstrap distribution. These bounds are not easily comparable to the stochastic bounds from Monte Carlo bootstrap confidence intervals. However, the narrowing of the mathematical bounds occurs at a faster rate than that of the Monte Carlo confidence intervals. This is demonstrated in Figure 3, which shows that for wide intervals, the Monte Carlo and discrete convolutional methods are comparable in terms of computation time. However, when more computation time is given, the discrete convolutional method is able to produce narrower intervals than the Monte Carlo bootstrap. Since the figure displays times and interval widths on the log scale, the differences are more obvious. Figure 3: A comparison of the two methods shows that the increase in computation time for smaller interval widths is much more dramatic with the Monte Carlo methods.

Monte Carlo confidence intervals narrow slowly as the computation time is increased. For the convolutional method the narrowing is approximately linear, in that doubling the computation time shrinks the interval by about one half. The Monte Carlo algorithm is often simple to implement with fast approximate answers, but when higher precision is desired, it slows considerably.

Figure 3 displays the computational time to approximate the bootstrap distribution of the sample mean for 20 observations. Since there is some natural variation in computation times this simulation was conducted 750 times and averaged together. The times for the Monte Carlo intervals are given by the solid curve and the times for the convolutional intervals are given by the dashed curve. The curves provide an approximation of the average computation times to achieve a given interval width. Notice that for the wider intervals the difference in computation time between the two methods is much less pronounced. In fact, when a wide interval is permissible, it is more computationally efficient to use the Monte Carlo bootstrap. However, the Monte Carlo method will only ever be able to generate stochastic bounds, whereas the discrete convolutional method is able to mathematically bound the estimate.

Referring again to Figure 3, the convolutional method’s computation times generally increase with the number of grid points selected along the support. Similarly, the Monte Carlo bootstrap also has narrower intervals as the number of samples (i.e., computation time) increases. However, the rate at which they narrow are clearly distinct, with the convolutional method converging faster.

Table 1 contains the discrete convolutional method times for the simulations in Figure 3. Note that the length of each support grid is a power of two, which takes advantage of this computational efficiency of the FFT. For samples that have high precision, a finer grid is needed to better capture the exact behavior of the bootstrapped distribution. Even with very fine grids the convolutional method is quite fast.

We also compare the convolutional method with saddlepoint approximations. We do this in several examples in the following section. Generally, the results from the convolutional method seem to be more accurate (i.e., they agree more closely with the Monte Carlo results) than the results from saddlepoint approximations.

## 4 Applications

A number of data applications have been chosen for inclusion in this article. These examples demonstrate the wide applicability of the convolutional method for bootstrapping. In each of the examples, saddlepoint approximations serve as the primary competitor to the convolutional method, and the Monte Carlo approach (with a very large number of samples) is considered to be the standard.

We argue that the convolutional approach has some distinct advantages over the saddlepoint method. First, the saddlepoint method is known to be quite accurate in the tails of distributions, but can be arbitrarily bad at specific points (Collins, 2009). The advantage of the convolutional method is that it can bound the error of the CDF at any point on the support. Another advantage is that the convolutional approach uses discrete distributions with compact support—which matches most bootstrap distributions. In contrast, the saddlepoint approximations are continuous and have infinite support. The final advantage is that the convolutional method is faster and easier to set up compared with the saddlepoint method, which requires at least calculating a second derivative of the cumulant generating function of the bootstrap distribution. We demonstrate these advantages in the following examples.

### 4.1 Example 2 from Davison and Hinkley (1988)

One example that lends itself well to comparison comes from Davison and Hinkley (1988). In this example, the authors look at the bootstrap distribution of the sample mean for the following 10 numbers:

Note that these values differ from those reported in the original paper. In the example from Davison and Hinkley (1988), the non-centered data were reported but were centered before performing the bootstrap; the centered observations are what have been reported above.

In the source article, 50,000 Monte Carlo samples were taken to approximate the distribution. With advances in computational capability, more precise approximations have been generated using 200M (million) Monte Carlo samples for this paper. We report the results from the 200M samples. Although Monte Carlo methods have no trouble estimating this distribution, to obtain this precision, it took 20.02 minutes of computation time.

The results using convolutional methods are very similar to the Monte Carlo estimates. However, the time needed for a machine to generate the convolutional results is practically instantaneous at just 0.0673 seconds. Table 2 shows the resulting mean absolute error (MAE) and mean squared error (MSE) between the two methods, which vanishes as the number of Monte Carlo samples increases. This implies the Monte Carlo results are converging to the convolutional results.

Naturally, the time needed to generate Monte Carlo samples is prone to change with simulation size. Nevertheless, we emphasize that the results from the convolutional method are exact.

The saddlepoint approximation suggested by Davison and Hinkley can be seen in Table 3. Compared to the Monte Carlo estimates the convolutional approach has an MAE of 0 and the saddlepoint’s MAE is 0.004. In general, these approximations are extremely close to the convolutional and Monte Carlo results. Though the convolutional and Monte Carlo results agree more often as indicated by the MAE.

Again, the results from the convolutional method described in this application are not an approximation. The distribution of is discrete and with no need to introduce discretization error, the quantiles described above are exact. Code for this and the following examples is found in the appendix.

### 4.2 Example 3 from Davison and Hinkley (1988)

Davison and Hinkley provide another interesting example in the same paper. This application studies the following twelve matched pairs differences with

While a bootstrap distribution of a statistic is still of interest, the convolution in question is somewhat different. Every value will be represented in the mean, however the sign has the ability to change. In other words, we are looking at the following convolution:

 ¯Y=Y1+Y2+...+Y1212

Where can take the values -4.5 and 4.5, can take the values -34.2 and 34.2, and so on. The probability of the positive or negative value appearing in a Monte Carlo sample is 0.5 for each of the random variables.

In our comparisons, we found the distribution for using six different Monte Carlo sample sizes. For each Monte Carlo sample, the 12

1 “sign” vector of -1 and 1’s was generated, -1 and 1 having the same probability of occurring in each element. The bootstrap sample was calculated as the mean of the element-wise product of the data and sign vectors. The time needed to find 1B (billion) Monte Carlo samples was more than 1.5 hours. As with the first example, the Monte Carlo estimates are very accurate. However, in much less time, exact quantiles were found using the convolutional method. Table

4 shows that as the number of Monte Carlo samples increases, the approximate Monte Carlo bootstrap distribution converges to the results of the convolutional method.

Table 5 compares the results of the three methods. The results were probabilities from the CDF at values specified in Davison and Hinkley (1988). The results from the convolutional method are exact, however, the values presented in Table 5 have been rounded to 5 decimal places. Measuring the distances from the Monte Carlo results we get an MAE of 0.00001 and 0.00364 for the convolutional and saddlepoint methods, respectively.

In this example, the convolutional method has better agreement with the Monte Carlo estimates than the saddlepoint method. There are a few quantiles whereon both the convolutional and saddlepoint methods agree, however the accuracy of the saddlepoint never exceeds the accuracy of the convolutional method. Once again, the computation time for the convolutional method is quite fast at 0.1130 seconds.

Generating the results from the convolutional method follows the procedure provided in Section 2 with one major exception. Rather than describing just one random variable, a random variable was needed for each of the pairs. From that point on, the basic theory of using multiplication in the Fourier domain to compute a convolution was utilized as in the previous example.

### 4.3 Laser Data Analysis from Balakrishnan and Qin (2019)

Balakrishnan and Qin (2019) provide an application using laser degradation data; this analysis lends itself well to comparison. Their application looks at cumulative laser device degradation using several different values as failure thresholds. In their analysis, they use a saddlepoint method to approximate the distribution of the 90th percentile of the time it takes to attain a cumulative degradation threshold. Here we use the phrase “failure time” to denote the random time until the designated failure threshold is achieved. Note that in the data, there are 16 degradation measurements per laser for 15 distinct lasers, resulting in a total of 240 measurements. The data can be found in Table C.17 of Meeker and Escobar (1998).

In order to find the 90th percentile of failure times, the following Monte Carlo resampling method was used. One of the 15 lasers was randomly chosen. Degradation measurements on the selected laser were randomly sampled (with replacement) until the cumulative degradation surpassed the threshold. Since the lasers were checked in increments of 250 hours, the exact failure times cannot typically be obtained. The final bootstrapped time was determined using linear interpolation on the last degradation measurement and the 250 hours. This process was repeated 40,000 times so that 40,000 bootstrapped failure times were collected. The 90th percentile of the bootstrapped times was recorded. This was then repeated 30,000 times, resulting in an approximation for the bootstrap distribution of the 90th percentile of failure times.

The convolutional method for this analysis also must take into account the specific laser effect, thus we extend the notation from Section 2.5 to accommodate distinct laser devices. We denote the measurement from the laser as (for and . These observations in turn define 15 discrete probability measures, with random variables , one for each laser device.

To obtain the laser failure time, , the method is the same as in Section 2.5. But accounting for all 15 lasers we must mix these CDFs. Thus for the bootstrap failure random variable, , its CDF is

 P(Z≤t)=11515∑i=1P(Zi≤t).

The distribution of combines the information of all lasers, but does not allow the data from one laser to affect the time of another. To find the 90th percentile of , a range of times was investigated.

The results of this analysis are included in Table 6. Notice that for each of the selected thresholds, the convolutional estimate is included in the Monte Carlo bootstrap interval. It is also interesting to note that the accuracy of the convolutional method seems to improve as the threshold increases. For a threshold of 1, these Monte Carlo estimates were found in 15.6 hours, however, 4.73 days were needed to find the estimates associated with a threshold of 10. The time required to find all ten estimates (upper and lower bounds) using the convolutional method was less than 6 minutes.

In this more complicated example, the convolutional method tends to agree with the Monte Carlo samples more often than the saddlepoint method. When treating the mean Monte Carlo results as the truth, the mean absolute error for the convolutional method is 0.038, while the MAE for the saddlepoint method is 8.935. It appears that the convolutional method has a clear advantage in terms of accuracy.

### 4.4 Numerical Example from Palayangoda et al. (2020)

In a more recent article, these laser data were analyzed again. The analysis was similar to that in Balakrishnan and Qin (2019), however in Palayangoda et al. (2020), data from individual lasers were treated as indistinguishable. More specifically, rather than only sampling from one laser at a time, sampling could occur from any laser degradation measurement.

Formally, the setup for the convolutional method in this example is shown in Section 2.5. Here all the data, for and , define the discrete random measure for the random variable .

In this example, the Monte Carlo estimates were calculated using 30,000 bootstrapped percentiles that were each calculated from 150,000 bootstrapped failure times. Table 7 shows the results from the three competing methods. When assuming the mean Monte Carlo estimates are the truth, the MAE for the convolutional method is 0.041 while the MAE for the saddlepoint method is 2.62. As in the previous comparison, the convolutional method appears to have an accuracy advantage over the saddlepoint method presented in Palayangoda et al. (2020).

In this simplified analysis, the Monte Carlo bootstrapping method had shorter computation times. It took 3.77 days to find Monte Carlo estimates for the distribution associated with a threshold of 10. Once again, the computation time required for the convolutional method is substantially shorter than the time needed for the Monte Carlo method, clocking in at just 16 seconds for all 10 thresholds and for both the lower and upper bounds. Although no computation times were reported for the saddlepoint method in Palayangoda et al. (2020), we assume they were much faster than the Monte Carlo method.

### 4.5 Asthma Application from Warr and Woodfield (2020)

One final example to consider is found in Warr and Woodfield (2020). The data for this application is found in the SemiMarkov package in R. The original study was conducted by observing asthma patients in France. The goal of this analysis is to approximate the first passage distribution from stage 1 to stage 3 asthma. See Figure 4 (adapted from Warr and Woodfield (2020)) for a graphical depiction of the states and transitions of asthma patients. Figure 4: Here the states and transitions of the asthma data are depicted. Each transition defines a unique distribution and occurs with fixed probability conditioned on the starting state. The figure has been adapted from Warr and Woodfield (2020).

The bootstrapping technique for this example is fairly straightforward. Each bootstrapped patient is assumed to start in stage 1. The next stage was chosen randomly with probabilities determined by sample proportions. A time was then randomly sampled from the observed times for that particular transition. This process was repeated until the bootstrapped patient transitioned into stage 3. In order to produce a comparison, one million Monte Carlo bootstrap samples were generated.

Furthermore, we computed two saddlepoint approximations for comparison. These approximations were also computationally quite fast, however a large degree of smoothing occurred. In order to get more accurate approximations, additional derivatives would need to be computed, making these approximations more difficult to obtain as accuracy increases.

The first of these two approximations is a second-order saddlepoint method (referred to as Saddlepoint 1). This method uses the approximation developed in Daniels (1954). The result of this method is an approximate pdf and so numerical integration was needed to convert to the CDF. The second saddlepoint method approximated the CDF directly but required an additional derivative to be computed. Saddlepoint 2 uses the approximation developed in Lugannani and Rice (1980).

The process of differentiation was made easier through the use of automatic differentiation. Specifically, the R package, numDeriv, was employed (Gilbert and Varadhan, 2019). This package computes accurate first and second order numerical derivatives. Since much of the difficulty in using saddlepoint approximations comes from the need to differentiate, automatic differentiation seems like an ideal tool to increase the feasibility of saddlepoint approximations using higher order terms to obtain more accuracy.

The analysis conducted in Warr and Woodfield differs in that the final result is a posterior curve for the CDF of first passage times. For this reason, the results here will not be directly compared to the results from that paper. However, we note that the methodology introduced in this article lends itself well to this application. In this case, the final result is an approximation of the CDF for the bootstrap distribution of first passage times. In this application the support of the bootstrap first passage distribution is unbounded, which complicates the mathematical bounding of the estimates. However, our approach still provides accurate estimates without appealing to more complex approaches.

Let denote the DFT sequence for the random variable defined by the time needed to transition from state to state , and denote the DFT sequence for the random variable defined by the first passage time needed to transition from state to state . Then the DFT for the first passage time from state 1 to state 3 is defined as:

 ~g1→3,k=p1~f1→3,k+(1−p1)~f1→2,k((1−p2)~f2→3,k+p1p2~f1→3,k~f2→1,k)1−(1−p1)p2~f1→2,k~f2→1,k, (4)

equivalent to Equation (5) in Warr and Woodfield. The theory for combining the DFTs in this way can be found in Pyke (1961). Here is the probability of transitioning directly to state 3 (instead of state 2), given the current state is state 1. is the probability of transitioning directly to state 1 (instead of state 3), given the current state is state 2.

Figure 5 displays the bootstrapped first-passage time distribution obtained through the discrete convolutional method. This distribution is compared to the results from the Monte Carlo and saddlepoint methods. Note, from a macro perspective, there is general agreement among the four methods. However, upon closer inspection it is clear the saddlepoint methods were not able to capture the discrete behavior inherent in this problem. The MAE between the Monte Carlo and convolutional methods was 0.001, while the MAE between the Monte Carlo and first and second saddlepoint methods was 0.0322 and 0.0216, respectively. Also note that less than 1 second was needed to compute the convolutional result, in contrast to the 11.3 seconds for the first saddlepoint result and 2.54 seconds for the second saddlepoint result. Approximately 6.53 hours of computation time were needed to generate the Monte Carlo samples. In producing the Monte Carlo samples, 100,000 first passage times were produced and then the quantiles of interest were recorded. This was iterated 15,000 times so that the uncertainty could be quantified. Figure 5: As seen from the image on the left, the four methods are in fairly close agreement. However, inspecting a close-up view of the approximations in the right image shows that the convolutional and Monte Carlo approaches are able to capture the same local behavior, but the two saddlepoint approaches impose a smoothed continuous approximation to the discrete distribution.

As a point of comparison, a few common quantiles were estimated from the first passage distribution using each method. The results are displayed in Table 8. Again, we point out that the Monte Carlo and convolutional methods agree quite closely.

For this application, the CDF was approximately bounded using the theory presented earlier. Table 9 presents the information for bounds generated on different grids. In each case, the grid started at 0 and ended at 30. Even for the grid of size the mean width of the error interval was less than 0.0001. As the grid becomes finer, the average difference in the error bounds appears to approach 0, as seen in Table 9. We note that due to the division in Equation 4 the bounds are not guaranteed, although we expect them to be quite accurate.

## 5 Conclusion

The convolutional bootstrap method provides an attractive alternative to the traditional Monte Carlo and saddlepoint methods. Due to the finite nature of samples, bootstrap distributions are discrete. Relying on this property, it is possible to use discrete computational algorithms to calculate the bootstrap distribution of some statistics. One simple application of the discrete convolutional method is in bootstrapping the sample mean. More interesting applications were demonstrated by bootstrapping failure times using degradation data and quantiles of first passage times.

Additionally, the discrete convolutional method makes it possible to find exact, or at least bounds on quantities of bootstrap distributions. The existence of mathematical error bounds provides a strong case for using the convolutional method. Even though Monte Carlo bootstrap methods provide confidence intervals, the intervals are stochastic and narrow at a slow rate.

The convolutional bootstrap method is very fast, due to the efficiency of the fast Fourier transform algorithm. We suppose our method is faster than the saddlepoint approach, and it is much faster than the Monte Carlo approach, for the achieved accuracy.

In this work we also point out that automatic differentiation could be a powerful tool when using saddlepoint approximations. Using this emerging numerical tool could allow saddlepoint approximations to be much faster to setup and easier to implement higher order approximations.

As a final note, the convolutional method is limited to statistics that are defined by sums of independent random variables. For example, we were unable to implement our method for the sample variance, since we know of no way to define it as a sum independent random variables. However, we demonstrated the convolutional method to be effective in bootstrapping means, first-passage times, and failure times from degradation data. The convolutional bootstrap method provides high levels of accuracy with relatively short computation times.

### Acknowledgements

The authors thank Dave Collins and Lynsie Warr for their content and editorial suggestions.

### Conflict of Interest

The authors have no conflicts of interest to declare that are relevant to the content of this article.

## References

• (1)
• Balakrishnan and Qin (2019) Balakrishnan, N. and Qin, C. (2019). Nonparametric evaluation of the first passage time of degradation processes, Applied Stochastic Models in Business and Industry 35(3): 571–590.
• Beyene (2001) Beyene, J. (2001). Uses of the fast Fourier transform (FFT) in exact statistical inference, PhD thesis, National Library of Canada.
• Butler and Bronson (2002) Butler, R. W. and Bronson, D. A. (2002). Bootstrapping survival times in stochastic systems by using saddlepoint approximations, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64(1): 31–49.
• Collins (2009) Collins, D. (2009). Nonparametric estimation of first passage time distributions in flowgraph models, PhD thesis, University of New Mexico.
• Cooley and Tukey (1965) Cooley, J. W. and Tukey, J. W. (1965). An algorithm for the machine calculation of complex fourier series, Mathematics of computation 19(90): 297–301.
• Daniels (1954) Daniels, H. E. (1954). Saddlepoint approximations in statistics, The Annals of Mathematical Statistics 25(4): 631–650.
• Davison and Hinkley (1988) Davison, A. C. and Hinkley, D. V. (1988). Saddlepoint approximations in resampling methods, Biometrika 75(3): 417–431.
• Efron (1979) Efron, B. (1979). Bootstrap methods: Another look at the jackknife, The Annals of Statistics 7(1): 1–26.
• Efron and Narasimhan (2020) Efron, B. and Narasimhan, B. (2020). The automatic construction of bootstrap confidence intervals, Journal of Computational and Graphical Statistics 29(3): 608–619.
• Efron and Tibshirani (1994) Efron, B. and Tibshirani, R. J. (1994). An introduction to the bootstrap, CRC press.
• Gilbert and Varadhan (2019) Gilbert, P. and Varadhan, R. (2019). numDeriv: Accurate Numerical Derivatives. R package version 2016.8-1.1.
https://CRAN.R-project.org/package=numDeriv
• Kern et al. (2003) Kern, J. W., McDonald, T. L., Amstrup, S. C., Durner, G. M. and Erickson, W. P. (2003). Using the bootstrap and fast fourier transform to estimate confidence intervals of 2d kernel densities, Environmental and Ecological Statistics 10(4): 405–418.
• Knight (1989) Knight, K. (1989). On the bootstrap of the sample mean in the infinite variance case, The Annals of Statistics 17(3): 1168–1175.
• Lugannani and Rice (1980) Lugannani, R. and Rice, S. (1980). Saddle point approximation for the distribution of the sum of independent random variables, Advances in applied probability 12(2): 475–490.
• Meeker and Escobar (1998) Meeker, W. Q. and Escobar, L. A. (1998). Statistical methods for reliability data, John Wiley & Sons.
• Palayangoda et al. (2020) Palayangoda, L. K., Ng, H. K. T. and Butler, R. W. (2020). Improved techniques for parametric and nonparametric evaluations of the first-passage time for degradation processes, Applied Stochastic Models in Business and Industry 36(4): 730–753.
• Pyke (1961) Pyke, R. (1961). Markov renewal processes with finitely many states, The Annals of Mathematical Statistics 32(4): 1243–1259.
• R Core Team (2021) R Core Team (2021). R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/
• Warr (2014) Warr, R. L. (2014). Numerical approximation of probability mass functions via the inverse discrete fourier transform, Methodology and Computing in Applied Probability 16(4): 1025–1038.
• Warr and Wight (2020) Warr, R. L. and Wight, C. J. (2020). Error bounds for cumulative distribution functions of convolutions via the discrete fourier transform, Methodology and Computing in Applied Probability 22(3): 881–904.
• Warr and Woodfield (2020) Warr, R. L. and Woodfield, T. B. (2020). Bayesian nonparametric estimation of first passage distributions in semi-markov processes, Applied Stochastic Models in Business and Industry 36(2): 237–250.

## Appendix A Code Appendix

### a.1 Example 2 from Davison and Hinkley (1988)

#The centered sample data
x <- c(-8.27,-7.46,-4.87,-2.87,-1.27,-0.67,-0.57,3.93,6.13,15.93)
#Transformed data
w <- round((x-min(x))/length(x),15)
#The support points
s <- (0:25000)/1000;
#Finding the pmf p on the support s
ind <- match(w,s)
p <- rep(0,length(s))
for (i in 1:length(x)) { p[ind[i]] <- p[ind[i]] + 1/length(x) }
#DFT of pmf defined by w
fp <- fft(p)
#DFT of ybar - min(x)
fm <- fp^length(x)
#pmf of ybar - min(x)
ybar <- Re(fft(fm,inverse=T))/length(s)
#Probabilities
probs <- c(0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.2,0.8,
0.9,0.95,0.99,0.995, 0.999, 0.9995, 0.9999)
quants <- numeric(length(probs))
for(i in 1:length(probs)){
quants[i] <- s[min(which(cumsum(ybar)>probs[i]))]+min(x)
}
quants #Results for convolutional method


### a.2 Example 3 from Davison and Hinkley (1988)

See Appendix B for a step-by-step description of the method and code in this application.

### a.3 Laser Data Analysis from Balakrishnan and Qin (2019)

#Find the data in Meeker and Escobar.
#The data are formatted such that it has 240 rows with the
# first 16 rows belong to laser-1 etc...
x <- round(read.csv(’laserdat.csv’)$increase,4) #Defining constants and the support N <- 2^18; d <- 250; s <- (0:(N-1))/10000 #Function to calculate CDF of Z using convolutional method cdf.Z <- function(Z,T){ m1 <- pmf.y(Z) CDF <- function(T) { Km1 <- floor(Z/d); index <- max(which(s <= T)); mix <- 0 for (i in 1:15) { fft.Zi <- fft(m1[i+15,])*fft(m1[i,])^Km1 mix <- mix + cumsum(Re(fft(fft.Zi,inverse=TRUE))/N)[index] } 1-mix/15 } CDF(T) } #Function to calculate pmf of y for some value of Z. # Each row of the output matrix is the pmf of y for the i-th laser # the last 15 rows includes a linear interpolation of the same # pmfs for the degradation in an interval less than 250 hours pmf.y <- function(Z) { pmf.yi <- matrix(0, nrow=30, ncol=N) for (i in 1:15) { pmf.yi[c(i,i+15),] <- yi.indices(i,Z) } pmf.yi } #Function which returns the indices (of the support) for the # point masses of y for the L-th laser. The first row of the # matrix is for degradations measured in 250 hour increments. # The 2nd row is a linear interpolation for smaller increments. # The last argument determines a lower or upper bound yi.indices <- function(L,Z){ prob250 <- numeric(N); probLI <- numeric(N); i <- 1:16 ind250 <- ceiling(x[i+16*(L-1)]*10000)+1 num250 <- table(ind250); prob250[unique(ind250)] <- num250/16 if (lower.bound) { indLI <- ceiling(x[i+16*(L-1)]*10000*(Z/d-floor(Z/d)))+1 } else { indLI <- floor(x[i+16*(L-1)]*10000*(Z/d-floor(Z/d)))+1 } numLI <- table(indLI); probLI[unique(indLI)] <- numLI/16 rbind(prob250, probLI) } #Defining a function such that the 90th percentile will be the root percent90 <- function(z,T) { cdf.Z(z,T)-0.9 } #0.9 quantiles for thresholds of 1,2,...,10 lower.quants <- rep(NA,10); upper.quants <- rep(NA,10); for (i in 1:10) { lower.bound=TRUE lower.quants[i] <- uniroot(percent90,c(250,6500),T=i)$root
lower.bound=FALSE
upper.quants[i] <- uniroot(percent90,c(250,6500),T=i)$root } round((lower.quants+upper.quants)/2,2)  ### a.4 Numerical Example from Palayangoda et al. (2020) # Run the code from the previous example first! #Function to calculate CDF of Z using convolutional method cdf.Z <- function(Z,T){ m1 <- pmf.y(Z) CDF <- function(T) { Km1 <- floor(Z/d); index <- max(which(s <= T)) fft.Z <- fft(m1[2,])*fft(m1[1,])^Km1 1-cumsum(Re(fft(fft.Z,inverse=TRUE))/N)[index] } CDF(T) } #Function to calculate pmf of y for some value of Z. # The 1st row of the output matrix is the pmf of y # The 2nd row includes a linear interpolation of the same # pmf for the degradation in an interval less than 250 hours pmf.y <- function(Z){ prob250 <- numeric(2^18); probLI <- numeric(2^18) for(laser in 1:15){ for(i in 1:16){ ind250 <-ceiling(x[i+16*(laser-1)]*10000)+1 prob250[ind250] <- prob250[ind250] + 1/16 if (lower.bound) { indLI <-ceiling(x[i+16*(laser-1)]*10000*(Z/d-floor(Z/d)))+1 } else { indLI <-floor(x[i+16*(laser-1)]*10000*(Z/d-floor(Z/d)))+1 } probLI[indLI] <- probLI[indLI] + 1/16 } } rbind(prob250, probLI)/15 } #0.9 quantiles for thresholds of 1,2,...,10 for (i in 1:10) { lower.bound=TRUE lower.quants[i] <- uniroot(percent90,c(250,6500),T=i)$root
lower.bound=FALSE
upper.quants[i] <- uniroot(percent90,c(250,6500),T=i)$root } round((lower.quants+upper.quants)/2,2)  ### a.5 Asthma Application from Warr and Woodfield (2020) #The data is contained in the SemiMarkov library library(SemiMarkov); data(asthma) #Initializing vectors and constants t12 <- numeric(); t13 <- numeric(); t21 <- numeric(); t23 <- numeric(); j <- 0; k <- 0; m <- 0; n <- 0 #Creating separate vectors for each of the transitions for(i in 1:nrow(asthma)){ if(asthma$state.h[i] == 1){
if(asthma$state.j[i] == 2){ j <- j+1; t12[j] <- asthma$time[i] }
if(asthma$state.j[i] == 3){ k <- k+1; t13[k] <- asthma$time[i] }
}
if(asthma$state.h[i] == 2){ if(asthma$state.j[i] == 1){
m <- m+1; t21[m] <- asthma$time[i] } if(asthma$state.j[i] == 3){
n <- n+1; t23[n] <- asthma\$time[i] }
} }
#The sample proportions serve as probability estimates
p12hat <- j/(j+k); p13hat <- k/(j+k)
p21hat <- m/(m+n); p23hat <- n/(m+n)
#An upper limit of 20 should capture virtually all the probability
s <- seq(0, 30, length.out=2^15)
#Defining the vectors of probabilities
pmf.yij <- function(times, s){
probs <- rep(0,length(s)); indices <- rep(0,length(times))
for(i in 1:length(times)){
indices[i] <- max(which(s<=times[i])) + upper.bound }
for (i in 1:length(indices)) {
probs[indices[i]] <- probs[indices[i]] + 1/length(times) }
probs
}
upper.bound=TRUE
pmf.12 <- pmf.yij(t12, s); pmf.13 <- pmf.yij(t13, s)
pmf.21 <- pmf.yij(t21, s); pmf.23 <- pmf.yij(t23, s)
#Employing the fft for each transition
ft.pmf.12 <- fft(p12hat*pmf.12); ft.pmf.13 <- fft(p13hat*pmf.13)
ft.pmf.21 <- fft(p21hat*pmf.21); ft.pmf.23 <- fft(p23hat*pmf.23)
#Convolution defined by transition from state 1 -> 3
fft.g13 <- ft.pmf.13+(ft.pmf.12*(ft.pmf.23+ft.pmf.13*ft.pmf.21)/
(1-ft.pmf.12*ft.pmf.21))
cdf.g13 <- cumsum(Re(fft(fft.g13,inverse=T))/length(s))
plot(s,cdf.g13,type="l") #Plot of the cdf
#Finding quantiles
qCon <- function(prob) s[max(which(cdf.g13<=prob))]
upper.quants <- c(qCon(.1),qCon(.25),qCon(.5),qCon(.75),qCon(.9))
#Repeating for the lower bounds
upper.bound=FALSE
pmf.12 <- pmf.yij(t12, s); pmf.13 <- pmf.yij(t13, s)
pmf.21 <- pmf.yij(t21, s); pmf.23 <- pmf.yij(t23, s)
ft.pmf.12 <- fft(p12hat*pmf.12); ft.pmf.13 <- fft(p13hat*pmf.13)
ft.pmf.21 <- fft(p21hat*pmf.21); ft.pmf.23 <- fft(p23hat*pmf.23)
fft.g13 <- ft.pmf.13+(ft.pmf.12*(ft.pmf.23+ft.pmf.13*ft.pmf.21)/
(1-ft.pmf.12*ft.pmf.21))
cdf.g13 <- cumsum(Re(fft(fft.g13,inverse=T))/length(s))
#Finding quantiles
lower.quants <- c(qCon(.1),qCon(.25),qCon(.5),qCon(.75),qCon(.9))
(upper.quants+lower.quants)/2


## Appendix B Example 3 from Davison and Hinkley (1988): Verbose Instructions

This section is meant to provide a step-by-step example of problem-solving with the convolutional method. The outline assumes that R is being used; however, this outline could be adjusted and applied to other programming languages.

• Start by storing the raw data in a vector.

x <- c(4.5,-34.2,7.4,12.6,-2.5,1.7,-34.0,7.3,15.4,-3.8,2.9,-4.2)


• Identify the convolution of interest. The convolution of interest for this application is the mean of 12 random variables, each with point mass of probability 0.5 at and , for .

• Create a separate vector for each of the random variables to be convolved. Each vector should contain the possible realizations for that RV. In this application, the vectors were started all at once by combining the data vectors, x, and -x to create a matrix with two columns.

n <- length(x); y <- cbind(x, -x);


• Create a vector that contains the support grid, s. Here the support was shifted to the right since fully positive supports are generally easier to work with. The support from 0 to 70 will allow for probability between -35 and 35.

s <- round(seq(0, 70, by=1/120),14)


• At this point multiple empty vectors were created that will be used in the following steps.

#Defining some data structures
p <- matrix(0, 12, length(s)); ind <- matrix(NA, 12, 2)
ind2 <- rep(NA,10); fft.p <- matrix(NA, 12, length(s))


• Before moving on, we want to make sure that the RVs described by the matrix y are the random variables that will be convolved. Here we shift the data so that it matches the support and then we divide by n. We also round at 14 decimal places to make sure the support numerically matches the adjusted data.

#The shifted and scaled data vectors
w <- round((y+35)/12,14)


• For each random variable described by the matrix w, we need to find which support points should have positive probability. In this step we are defining the pmf for each random variable and inputting it as a row in the matrix p.

#Defining the pmf of each w on s
for (i in 1:n) { ind[i,] <- match(w[i,],s) }
for (j in 1:n) { for (i in 1:2){
p[j,ind[j,i]] <- p[j,ind[j,i]] + 1/2
} }


• Now the FFT is applied to each random variable by using the fft function on each row of the matrix, p.

#Computing the fft for each pmf
for (i in 1:nrow(p)) { fft.p[i,] <- fft(p[i,]) }


• At this point, we perform the convolution of interest using the multiplicative property of the DFT.

#Finding the fft for the pmf of ybar
fft.ybar <- apply(fft.p,2,prod)


• Finally, we are able to invert the FFT to find the pmf for the convolution of interest. Notice that in R we must divide by the size of the support in order to normalize the pmf. We are also finding the cumulative sum in order to obtain the CDF. The quantiles were obtained in order to compare to the results from the source article.

#The CDF of ybar
cdf.ybar <- cumsum(Re(fft(fft.ybar,inverse=T))/length(s))
#Finding the CDF values at specific support points
quants <- c(-10.77,-10.32,-8.97,-8.53,-7.63,-6.28,-4.04,-2.24,-.9,0)
for (i in 1:10) { ind2[i] <- max(which(s-35 <= quants[i])) }
cdf.ybar[ind2]