The Direct Convolution (DC) approach is requested with
method = "Convolve".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
dpbinom(NULL, pp, wt, "Convolve")
#> [1] 3.574462e-35 1.120280e-32 1.685184e-30 1.620524e-28 1.119523e-26
#> [6] 5.920060e-25 2.493263e-23 8.591850e-22 2.470125e-20 6.011429e-19
#> [11] 1.252345e-17 2.253115e-16 3.525477e-15 4.825171e-14 5.803728e-13
#> [16] 6.158735e-12 5.784692e-11 4.822437e-10 3.576566e-09 2.364563e-08
#> [21] 1.395965e-07 7.370448e-07 3.484836e-06 1.477208e-05 5.619632e-05
#> [26] 1.920240e-04 5.897928e-04 1.629272e-03 4.049768e-03 9.060183e-03
#> [31] 1.824629e-02 3.307754e-02 5.396724e-02 7.921491e-02 1.045505e-01
#> [36] 1.239854e-01 1.319896e-01 1.259938e-01 1.077029e-01 8.232174e-02
#> [41] 5.616422e-02 3.413623e-02 1.844304e-02 8.835890e-03 3.743554e-03
#> [46] 1.398320e-03 4.589049e-04 1.318064e-04 3.298425e-05 7.154649e-06
#> [51] 1.337083e-06 2.137543e-07 2.898296e-08 3.298587e-09 3.110922e-10
#> [56] 2.392070e-11 1.468267e-12 6.991155e-14 2.478218e-15 6.130807e-17
#> [61] 9.411166e-19 6.727527e-21
ppbinom(NULL, pp, wt, "Convolve")
#> [1] 3.574462e-35 1.123854e-32 1.696423e-30 1.637488e-28 1.135898e-26
#> [6] 6.033650e-25 2.553600e-23 8.847210e-22 2.558597e-20 6.267289e-19
#> [11] 1.315018e-17 2.384617e-16 3.763939e-15 5.201565e-14 6.323884e-13
#> [16] 6.791123e-12 6.463805e-11 5.468818e-10 4.123448e-09 2.776908e-08
#> [21] 1.673656e-07 9.044104e-07 4.389247e-06 1.916133e-05 7.535765e-05
#> [26] 2.673817e-04 8.571745e-04 2.486446e-03 6.536215e-03 1.559640e-02
#> [31] 3.384269e-02 6.692022e-02 1.208875e-01 2.001024e-01 3.046529e-01
#> [36] 4.286383e-01 5.606280e-01 6.866217e-01 7.943246e-01 8.766463e-01
#> [41] 9.328105e-01 9.669468e-01 9.853898e-01 9.942257e-01 9.979692e-01
#> [46] 9.993676e-01 9.998265e-01 9.999583e-01 9.999913e-01 9.999984e-01
#> [51] 9.999998e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [56] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [61] 1.000000e+00 1.000000e+00The Divide & Conquer FFT Tree Convolution (DC-FFT)
approach is requested with method = "DivideFFT".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
dpbinom(NULL, pp, wt, "DivideFFT")
#> [1] 3.574462e-35 1.120280e-32 1.685184e-30 1.620524e-28 1.119523e-26
#> [6] 5.920060e-25 2.493263e-23 8.591850e-22 2.470125e-20 6.011429e-19
#> [11] 1.252345e-17 2.253115e-16 3.525477e-15 4.825171e-14 5.803728e-13
#> [16] 6.158735e-12 5.784692e-11 4.822437e-10 3.576566e-09 2.364563e-08
#> [21] 1.395965e-07 7.370448e-07 3.484836e-06 1.477208e-05 5.619632e-05
#> [26] 1.920240e-04 5.897928e-04 1.629272e-03 4.049768e-03 9.060183e-03
#> [31] 1.824629e-02 3.307754e-02 5.396724e-02 7.921491e-02 1.045505e-01
#> [36] 1.239854e-01 1.319896e-01 1.259938e-01 1.077029e-01 8.232174e-02
#> [41] 5.616422e-02 3.413623e-02 1.844304e-02 8.835890e-03 3.743554e-03
#> [46] 1.398320e-03 4.589049e-04 1.318064e-04 3.298425e-05 7.154649e-06
#> [51] 1.337083e-06 2.137543e-07 2.898296e-08 3.298587e-09 3.110922e-10
#> [56] 2.392070e-11 1.468267e-12 6.991155e-14 2.478218e-15 6.130807e-17
#> [61] 9.411166e-19 6.727527e-21
ppbinom(NULL, pp, wt, "DivideFFT")
#> [1] 3.574462e-35 1.123854e-32 1.696423e-30 1.637488e-28 1.135898e-26
#> [6] 6.033650e-25 2.553600e-23 8.847210e-22 2.558597e-20 6.267289e-19
#> [11] 1.315018e-17 2.384617e-16 3.763939e-15 5.201565e-14 6.323884e-13
#> [16] 6.791123e-12 6.463805e-11 5.468818e-10 4.123448e-09 2.776908e-08
#> [21] 1.673656e-07 9.044104e-07 4.389247e-06 1.916133e-05 7.535765e-05
#> [26] 2.673817e-04 8.571745e-04 2.486446e-03 6.536215e-03 1.559640e-02
#> [31] 3.384269e-02 6.692022e-02 1.208875e-01 2.001024e-01 3.046529e-01
#> [36] 4.286383e-01 5.606280e-01 6.866217e-01 7.943246e-01 8.766463e-01
#> [41] 9.328105e-01 9.669468e-01 9.853898e-01 9.942257e-01 9.979692e-01
#> [46] 9.993676e-01 9.998265e-01 9.999583e-01 9.999913e-01 9.999984e-01
#> [51] 9.999998e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [56] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [61] 1.000000e+00 1.000000e+00By design, as proposed by Biscarri, Zhao & Brunner (2018), its results are identical to the DC procedure, if \(n \leq 750\). Thus, differences can be observed for larger \(n > 750\):
set.seed(1)
pp1 <- runif(751)
pp2 <- pp1[1:750]
sum(abs(dpbinom(NULL, pp2, method = "DivideFFT") - dpbinom(NULL, pp2, method = "Convolve")))
#> [1] 0
sum(abs(dpbinom(NULL, pp1, method = "DivideFFT") - dpbinom(NULL, pp1, method = "Convolve")))
#> [1] 0The reason is that the DC-FFT method splits the input
probs vector into as equally sized parts as possible and
computes their distributions separately with the DC approach. The
results of the portions are then convoluted by means of the Fast Fourier
Transformation. As proposed by Biscarri, Zhao &
Brunner (2018), no splitting is done for \(n \leq 750\). In addition, the DC-FFT
procedure does not produce probabilities \(\leq 5.55e\text{-}17\), i.e. smaller values
are rounded off to 0, if \(n >
750\), whereas the smallest possible result of the DC algorithm
is \(\sim 1e\text{-}323\). This is most
likely caused by the used FFTW3 library.
The Discrete Fourier Transformation of the Characteristic
Function (DFT-CF) approach is requested with
method = "Characteristic".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
dpbinom(NULL, pp, wt, "Characteristic")
#> [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [11] 0.000000e+00 2.238353e-16 3.549132e-15 4.829828e-14 5.804377e-13
#> [16] 6.158818e-12 5.784702e-11 4.822438e-10 3.576566e-09 2.364563e-08
#> [21] 1.395965e-07 7.370448e-07 3.484836e-06 1.477208e-05 5.619632e-05
#> [26] 1.920240e-04 5.897928e-04 1.629272e-03 4.049768e-03 9.060183e-03
#> [31] 1.824629e-02 3.307754e-02 5.396724e-02 7.921491e-02 1.045505e-01
#> [36] 1.239854e-01 1.319896e-01 1.259938e-01 1.077029e-01 8.232174e-02
#> [41] 5.616422e-02 3.413623e-02 1.844304e-02 8.835890e-03 3.743554e-03
#> [46] 1.398320e-03 4.589049e-04 1.318064e-04 3.298425e-05 7.154649e-06
#> [51] 1.337083e-06 2.137543e-07 2.898296e-08 3.298587e-09 3.110923e-10
#> [56] 2.392079e-11 1.468354e-12 6.994931e-14 2.513558e-15 0.000000e+00
#> [61] 0.000000e+00 0.000000e+00
ppbinom(NULL, pp, wt, "Characteristic")
#> [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [11] 0.000000e+00 2.238353e-16 3.772968e-15 5.207125e-14 6.325089e-13
#> [16] 6.791327e-12 6.463834e-11 5.468822e-10 4.123448e-09 2.776908e-08
#> [21] 1.673656e-07 9.044104e-07 4.389247e-06 1.916133e-05 7.535765e-05
#> [26] 2.673817e-04 8.571745e-04 2.486446e-03 6.536215e-03 1.559640e-02
#> [31] 3.384269e-02 6.692022e-02 1.208875e-01 2.001024e-01 3.046529e-01
#> [36] 4.286383e-01 5.606280e-01 6.866217e-01 7.943246e-01 8.766463e-01
#> [41] 9.328105e-01 9.669468e-01 9.853898e-01 9.942257e-01 9.979692e-01
#> [46] 9.993676e-01 9.998265e-01 9.999583e-01 9.999913e-01 9.999984e-01
#> [51] 9.999998e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [56] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [61] 1.000000e+00 1.000000e+00As can be seen, the DFT-CF procedure does not produce probabilities \(\leq 2.22e\text{-}16\), i.e. smaller values are rounded off to 0, most likely due to the used FFTW3 library.
The Recursive Formula (RF) approach is requested with
method = "Recursive".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
dpbinom(NULL, pp, wt, "Recursive")
#> [1] 3.574462e-35 1.120280e-32 1.685184e-30 1.620524e-28 1.119523e-26
#> [6] 5.920060e-25 2.493263e-23 8.591850e-22 2.470125e-20 6.011429e-19
#> [11] 1.252345e-17 2.253115e-16 3.525477e-15 4.825171e-14 5.803728e-13
#> [16] 6.158735e-12 5.784692e-11 4.822437e-10 3.576566e-09 2.364563e-08
#> [21] 1.395965e-07 7.370448e-07 3.484836e-06 1.477208e-05 5.619632e-05
#> [26] 1.920240e-04 5.897928e-04 1.629272e-03 4.049768e-03 9.060183e-03
#> [31] 1.824629e-02 3.307754e-02 5.396724e-02 7.921491e-02 1.045505e-01
#> [36] 1.239854e-01 1.319896e-01 1.259938e-01 1.077029e-01 8.232174e-02
#> [41] 5.616422e-02 3.413623e-02 1.844304e-02 8.835890e-03 3.743554e-03
#> [46] 1.398320e-03 4.589049e-04 1.318064e-04 3.298425e-05 7.154649e-06
#> [51] 1.337083e-06 2.137543e-07 2.898296e-08 3.298587e-09 3.110922e-10
#> [56] 2.392070e-11 1.468267e-12 6.991155e-14 2.478218e-15 6.130807e-17
#> [61] 9.411166e-19 6.727527e-21
ppbinom(NULL, pp, wt, "Recursive")
#> [1] 3.574462e-35 1.123854e-32 1.696423e-30 1.637488e-28 1.135898e-26
#> [6] 6.033650e-25 2.553600e-23 8.847210e-22 2.558597e-20 6.267289e-19
#> [11] 1.315018e-17 2.384617e-16 3.763939e-15 5.201565e-14 6.323884e-13
#> [16] 6.791123e-12 6.463805e-11 5.468818e-10 4.123448e-09 2.776908e-08
#> [21] 1.673656e-07 9.044104e-07 4.389247e-06 1.916133e-05 7.535765e-05
#> [26] 2.673817e-04 8.571745e-04 2.486446e-03 6.536215e-03 1.559640e-02
#> [31] 3.384269e-02 6.692022e-02 1.208875e-01 2.001024e-01 3.046529e-01
#> [36] 4.286383e-01 5.606280e-01 6.866217e-01 7.943246e-01 8.766463e-01
#> [41] 9.328105e-01 9.669468e-01 9.853898e-01 9.942257e-01 9.979692e-01
#> [46] 9.993676e-01 9.998265e-01 9.999583e-01 9.999913e-01 9.999984e-01
#> [51] 9.999998e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [56] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [61] 1.000000e+00 1.000000e+00Obviously, the RF procedure does produce probabilities \(\leq 5.55e\text{-}17\), because it does not rely on the FFTW3 library. Furthermore, it yields the same results as the DC method.
To assess the performance of the exact procedures, we use the
microbenchmark package. Each algorithm has to calculate the
PMF repeatedly based on random probability vectors. The run times are
then summarized in a table that presents, among other statistics, their
minima, maxima and means. The following results were recorded on an AMD
Ryzen 9 5900X with 64 GiB of RAM and Manjaro KDE Linux (26.0.2).
library(microbenchmark)
set.seed(1)
f1 <- function() dpbinom(NULL, runif(6000), method = "DivideFFT")
f2 <- function() dpbinom(NULL, runif(6000), method = "Convolve")
f3 <- function() dpbinom(NULL, runif(6000), method = "Recursive")
f4 <- function() dpbinom(NULL, runif(6000), method = "Characteristic")
microbenchmark(f1(), f2(), f3(), f4(), times = 51)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> f1() 19.71430 19.80681 20.04784 19.91380 20.00289 22.92915 51
#> f2() 49.22795 49.46415 49.68863 49.60254 49.83205 51.20837 51
#> f3() 123.76136 124.02034 124.31687 124.20418 124.42785 127.07334 51
#> f4() 140.34830 143.47662 144.40956 144.51428 145.16204 147.47353 51Clearly, the DC-FFT procedure is the fastest, followed by DC, RF and DFT-CF methods.
The Generalized Direct Convolution (G-DC) approach is
requested with method = "Convolve".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
va <- sample(0:10, 10, TRUE)
vb <- sample(0:10, 10, TRUE)
dgpbinom(NULL, pp, va, vb, wt, "Convolve")
#> [1] 1.140600e-31 5.349930e-30 1.164698e-28 1.572037e-27 1.491024e-26
#> [6] 1.077204e-25 6.336147e-25 3.215011e-24 1.466295e-23 6.127671e-23
#> [11] 2.363402e-22 8.484857e-22 2.866109e-21 9.171228e-21 2.788507e-20
#> [16] 8.091940e-20 2.254155e-19 6.051395e-19 1.570129e-18 3.953458e-18
#> [21] 9.696098e-18 2.321913e-17 5.442392e-17 1.251302e-16 2.824507e-16
#> [26] 6.264454e-16 1.366745e-15 2.934598e-15 6.203639e-15 1.292697e-14
#> [31] 2.657759e-14 5.394727e-14 1.081983e-13 2.144873e-13 4.201625e-13
#> [36] 8.135609e-13 1.557745e-12 2.949821e-12 5.527695e-12 1.025815e-11
#> [41] 1.885777e-11 3.434641e-11 6.196981e-11 1.106787e-10 1.956340e-10
#> [46] 3.425394e-10 5.948077e-10 1.025224e-09 1.753751e-09 2.972596e-09
#> [51] 4.985314e-09 8.275458e-09 1.362195e-08 2.227979e-08 3.622799e-08
#> [56] 5.845270e-08 9.332219e-08 1.473012e-07 2.302797e-07 3.576650e-07
#> [61] 5.529336e-07 8.496291e-07 1.292864e-06 1.943382e-06 2.888042e-06
#> [66] 4.257944e-06 6.248675e-06 9.128095e-06 1.322640e-05 1.893515e-05
#> [71] 2.675612e-05 3.741507e-05 5.199255e-05 7.194684e-05 9.895330e-05
#> [76] 1.347017e-04 1.809349e-04 2.399008e-04 3.150314e-04 4.112231e-04
#> [81] 5.341537e-04 6.888863e-04 8.788234e-04 1.106198e-03 1.374340e-03
#> [86] 1.690272e-03 2.065290e-03 2.511885e-03 3.037800e-03 3.641214e-03
#> [91] 4.311837e-03 5.039293e-03 5.824625e-03 6.686091e-03 7.651765e-03
#> [96] 8.740859e-03 9.945159e-03 1.122411e-02 1.252016e-02 1.378863e-02
#> [101] 1.502576e-02 1.627450e-02 1.759663e-02 1.902489e-02 2.052786e-02
#> [106] 2.201243e-02 2.336424e-02 2.450429e-02 2.543095e-02 2.622065e-02
#> [111] 2.697857e-02 2.776636e-02 2.855637e-02 2.924236e-02 2.969655e-02
#> [116] 2.983772e-02 2.967384e-02 2.929746e-02 2.883252e-02 2.836282e-02
#> [121] 2.788971e-02 2.734351e-02 2.663438e-02 2.570794e-02 2.457639e-02
#> [126] 2.331289e-02 2.201380e-02 2.075053e-02 1.954176e-02 1.836001e-02
#> [131] 1.716200e-02 1.592047e-02 1.464084e-02 1.335803e-02 1.211826e-02
#> [136] 1.095708e-02 9.886542e-03 8.897658e-03 7.972694e-03 7.098018e-03
#> [141] 6.270583e-03 5.496952e-03 4.787457e-03 4.149442e-03 3.583427e-03
#> [146] 3.083701e-03 2.641746e-03 2.249767e-03 1.902455e-03 1.596805e-03
#> [151] 1.330879e-03 1.102475e-03 9.084265e-04 7.447312e-04 6.071616e-04
#> [156] 4.918629e-04 3.956251e-04 3.158260e-04 2.502339e-04 1.968330e-04
#> [161] 1.537458e-04 1.192445e-04 9.179821e-05 7.010494e-05 5.308547e-05
#> [166] 3.984854e-05 2.965115e-05 2.187013e-05 1.598631e-05 1.157497e-05
#> [171] 8.295941e-06 5.881266e-06 4.121776e-06 2.854642e-06 1.953341e-06
#> [176] 1.320224e-06 8.809465e-07 5.799307e-07 3.763587e-07 2.406488e-07
#> [181] 1.515662e-07 9.401686e-08 5.742327e-08 3.451481e-08 2.039831e-08
#> [186] 1.184350e-08 6.751380e-09 3.777327e-09 2.073644e-09 1.116337e-09
#> [191] 5.887148e-10 3.036829e-10 1.529887e-10 7.516829e-11 3.598151e-11
#> [196] 1.676154e-11 7.585978e-12 3.326429e-12 1.407527e-12 5.717370e-13
#> [201] 2.216349e-13 8.149241e-14 2.824954e-14 9.179165e-15 2.780017e-15
#> [206] 7.803525e-16 2.018046e-16 4.775552e-17 1.025798e-17 1.979767e-18
#> [211] 3.386554e-19 5.038594e-20 6.336865e-21 6.424747e-22 4.821385e-23
#> [216] 2.108301e-24
pgpbinom(NULL, pp, va, vb, wt, "Convolve")
#> [1] 1.140600e-31 5.463990e-30 1.219337e-28 1.693971e-27 1.660421e-26
#> [6] 1.243246e-25 7.579393e-25 3.972950e-24 1.863590e-23 7.991261e-23
#> [11] 3.162528e-22 1.164739e-21 4.030847e-21 1.320208e-20 4.108715e-20
#> [16] 1.220065e-19 3.474220e-19 9.525615e-19 2.522691e-18 6.476149e-18
#> [21] 1.617225e-17 3.939138e-17 9.381530e-17 2.189455e-16 5.013962e-16
#> [26] 1.127842e-15 2.494586e-15 5.429184e-15 1.163282e-14 2.455979e-14
#> [31] 5.113739e-14 1.050847e-13 2.132829e-13 4.277703e-13 8.479327e-13
#> [36] 1.661494e-12 3.219239e-12 6.169059e-12 1.169675e-11 2.195491e-11
#> [41] 4.081268e-11 7.515909e-11 1.371289e-10 2.478076e-10 4.434415e-10
#> [46] 7.859810e-10 1.380789e-09 2.406013e-09 4.159763e-09 7.132360e-09
#> [51] 1.211767e-08 2.039313e-08 3.401508e-08 5.629487e-08 9.252285e-08
#> [56] 1.509756e-07 2.442977e-07 3.915989e-07 6.218786e-07 9.795436e-07
#> [61] 1.532477e-06 2.382106e-06 3.674970e-06 5.618352e-06 8.506394e-06
#> [66] 1.276434e-05 1.901301e-05 2.814111e-05 4.136751e-05 6.030266e-05
#> [71] 8.705877e-05 1.244738e-04 1.764664e-04 2.484132e-04 3.473665e-04
#> [76] 4.820683e-04 6.630032e-04 9.029039e-04 1.217935e-03 1.629158e-03
#> [81] 2.163312e-03 2.852198e-03 3.731022e-03 4.837220e-03 6.211560e-03
#> [86] 7.901832e-03 9.967122e-03 1.247901e-02 1.551681e-02 1.915802e-02
#> [91] 2.346986e-02 2.850915e-02 3.433378e-02 4.101987e-02 4.867163e-02
#> [96] 5.741249e-02 6.735765e-02 7.858176e-02 9.110192e-02 1.048906e-01
#> [101] 1.199163e-01 1.361908e-01 1.537874e-01 1.728123e-01 1.933402e-01
#> [106] 2.153526e-01 2.387169e-01 2.632211e-01 2.886521e-01 3.148727e-01
#> [111] 3.418513e-01 3.696177e-01 3.981740e-01 4.274164e-01 4.571130e-01
#> [116] 4.869507e-01 5.166245e-01 5.459220e-01 5.747545e-01 6.031173e-01
#> [121] 6.310070e-01 6.583505e-01 6.849849e-01 7.106929e-01 7.352692e-01
#> [126] 7.585821e-01 7.805959e-01 8.013465e-01 8.208882e-01 8.392482e-01
#> [131] 8.564102e-01 8.723307e-01 8.869715e-01 9.003296e-01 9.124478e-01
#> [136] 9.234049e-01 9.332914e-01 9.421891e-01 9.501618e-01 9.572598e-01
#> [141] 9.635304e-01 9.690273e-01 9.738148e-01 9.779642e-01 9.815477e-01
#> [146] 9.846314e-01 9.872731e-01 9.895229e-01 9.914253e-01 9.930221e-01
#> [151] 9.943530e-01 9.954555e-01 9.963639e-01 9.971087e-01 9.977158e-01
#> [156] 9.982077e-01 9.986033e-01 9.989191e-01 9.991694e-01 9.993662e-01
#> [161] 9.995199e-01 9.996392e-01 9.997310e-01 9.998011e-01 9.998542e-01
#> [166] 9.998940e-01 9.999237e-01 9.999455e-01 9.999615e-01 9.999731e-01
#> [171] 9.999814e-01 9.999873e-01 9.999914e-01 9.999943e-01 9.999962e-01
#> [176] 9.999975e-01 9.999984e-01 9.999990e-01 9.999994e-01 9.999996e-01
#> [181] 9.999998e-01 9.999999e-01 9.999999e-01 1.000000e+00 1.000000e+00
#> [186] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [191] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [196] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [201] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [206] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [211] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [216] 1.000000e+00The Generalized Divide & Conquer FFT Tree Convolution
(G-DC-FFT) approach is requested with
method = "DivideFFT".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
va <- sample(0:10, 10, TRUE)
vb <- sample(0:10, 10, TRUE)
dgpbinom(NULL, pp, va, vb, wt, "DivideFFT")
#> [1] 1.140600e-31 5.349930e-30 1.164698e-28 1.572037e-27 1.491024e-26
#> [6] 1.077204e-25 6.336147e-25 3.215011e-24 1.466295e-23 6.127671e-23
#> [11] 2.363402e-22 8.484857e-22 2.866109e-21 9.171228e-21 2.788507e-20
#> [16] 8.091940e-20 2.254155e-19 6.051395e-19 1.570129e-18 3.953458e-18
#> [21] 9.696098e-18 2.321913e-17 5.442392e-17 1.251302e-16 2.824507e-16
#> [26] 6.264454e-16 1.366745e-15 2.934598e-15 6.203639e-15 1.292697e-14
#> [31] 2.657759e-14 5.394727e-14 1.081983e-13 2.144873e-13 4.201625e-13
#> [36] 8.135609e-13 1.557745e-12 2.949821e-12 5.527695e-12 1.025815e-11
#> [41] 1.885777e-11 3.434641e-11 6.196981e-11 1.106787e-10 1.956340e-10
#> [46] 3.425394e-10 5.948077e-10 1.025224e-09 1.753751e-09 2.972596e-09
#> [51] 4.985314e-09 8.275458e-09 1.362195e-08 2.227979e-08 3.622799e-08
#> [56] 5.845270e-08 9.332219e-08 1.473012e-07 2.302797e-07 3.576650e-07
#> [61] 5.529336e-07 8.496291e-07 1.292864e-06 1.943382e-06 2.888042e-06
#> [66] 4.257944e-06 6.248675e-06 9.128095e-06 1.322640e-05 1.893515e-05
#> [71] 2.675612e-05 3.741507e-05 5.199255e-05 7.194684e-05 9.895330e-05
#> [76] 1.347017e-04 1.809349e-04 2.399008e-04 3.150314e-04 4.112231e-04
#> [81] 5.341537e-04 6.888863e-04 8.788234e-04 1.106198e-03 1.374340e-03
#> [86] 1.690272e-03 2.065290e-03 2.511885e-03 3.037800e-03 3.641214e-03
#> [91] 4.311837e-03 5.039293e-03 5.824625e-03 6.686091e-03 7.651765e-03
#> [96] 8.740859e-03 9.945159e-03 1.122411e-02 1.252016e-02 1.378863e-02
#> [101] 1.502576e-02 1.627450e-02 1.759663e-02 1.902489e-02 2.052786e-02
#> [106] 2.201243e-02 2.336424e-02 2.450429e-02 2.543095e-02 2.622065e-02
#> [111] 2.697857e-02 2.776636e-02 2.855637e-02 2.924236e-02 2.969655e-02
#> [116] 2.983772e-02 2.967384e-02 2.929746e-02 2.883252e-02 2.836282e-02
#> [121] 2.788971e-02 2.734351e-02 2.663438e-02 2.570794e-02 2.457639e-02
#> [126] 2.331289e-02 2.201380e-02 2.075053e-02 1.954176e-02 1.836001e-02
#> [131] 1.716200e-02 1.592047e-02 1.464084e-02 1.335803e-02 1.211826e-02
#> [136] 1.095708e-02 9.886542e-03 8.897658e-03 7.972694e-03 7.098018e-03
#> [141] 6.270583e-03 5.496952e-03 4.787457e-03 4.149442e-03 3.583427e-03
#> [146] 3.083701e-03 2.641746e-03 2.249767e-03 1.902455e-03 1.596805e-03
#> [151] 1.330879e-03 1.102475e-03 9.084265e-04 7.447312e-04 6.071616e-04
#> [156] 4.918629e-04 3.956251e-04 3.158260e-04 2.502339e-04 1.968330e-04
#> [161] 1.537458e-04 1.192445e-04 9.179821e-05 7.010494e-05 5.308547e-05
#> [166] 3.984854e-05 2.965115e-05 2.187013e-05 1.598631e-05 1.157497e-05
#> [171] 8.295941e-06 5.881266e-06 4.121776e-06 2.854642e-06 1.953341e-06
#> [176] 1.320224e-06 8.809465e-07 5.799307e-07 3.763587e-07 2.406488e-07
#> [181] 1.515662e-07 9.401686e-08 5.742327e-08 3.451481e-08 2.039831e-08
#> [186] 1.184350e-08 6.751380e-09 3.777327e-09 2.073644e-09 1.116337e-09
#> [191] 5.887148e-10 3.036829e-10 1.529887e-10 7.516829e-11 3.598151e-11
#> [196] 1.676154e-11 7.585978e-12 3.326429e-12 1.407527e-12 5.717370e-13
#> [201] 2.216349e-13 8.149241e-14 2.824954e-14 9.179165e-15 2.780017e-15
#> [206] 7.803525e-16 2.018046e-16 4.775552e-17 1.025798e-17 1.979767e-18
#> [211] 3.386554e-19 5.038594e-20 6.336865e-21 6.424747e-22 4.821385e-23
#> [216] 2.108301e-24
pgpbinom(NULL, pp, va, vb, wt, "DivideFFT")
#> [1] 1.140600e-31 5.463990e-30 1.219337e-28 1.693971e-27 1.660421e-26
#> [6] 1.243246e-25 7.579393e-25 3.972950e-24 1.863590e-23 7.991261e-23
#> [11] 3.162528e-22 1.164739e-21 4.030847e-21 1.320208e-20 4.108715e-20
#> [16] 1.220065e-19 3.474220e-19 9.525615e-19 2.522691e-18 6.476149e-18
#> [21] 1.617225e-17 3.939138e-17 9.381530e-17 2.189455e-16 5.013962e-16
#> [26] 1.127842e-15 2.494586e-15 5.429184e-15 1.163282e-14 2.455979e-14
#> [31] 5.113739e-14 1.050847e-13 2.132829e-13 4.277703e-13 8.479327e-13
#> [36] 1.661494e-12 3.219239e-12 6.169059e-12 1.169675e-11 2.195491e-11
#> [41] 4.081268e-11 7.515909e-11 1.371289e-10 2.478076e-10 4.434415e-10
#> [46] 7.859810e-10 1.380789e-09 2.406013e-09 4.159763e-09 7.132360e-09
#> [51] 1.211767e-08 2.039313e-08 3.401508e-08 5.629487e-08 9.252285e-08
#> [56] 1.509756e-07 2.442977e-07 3.915989e-07 6.218786e-07 9.795436e-07
#> [61] 1.532477e-06 2.382106e-06 3.674970e-06 5.618352e-06 8.506394e-06
#> [66] 1.276434e-05 1.901301e-05 2.814111e-05 4.136751e-05 6.030266e-05
#> [71] 8.705877e-05 1.244738e-04 1.764664e-04 2.484132e-04 3.473665e-04
#> [76] 4.820683e-04 6.630032e-04 9.029039e-04 1.217935e-03 1.629158e-03
#> [81] 2.163312e-03 2.852198e-03 3.731022e-03 4.837220e-03 6.211560e-03
#> [86] 7.901832e-03 9.967122e-03 1.247901e-02 1.551681e-02 1.915802e-02
#> [91] 2.346986e-02 2.850915e-02 3.433378e-02 4.101987e-02 4.867163e-02
#> [96] 5.741249e-02 6.735765e-02 7.858176e-02 9.110192e-02 1.048906e-01
#> [101] 1.199163e-01 1.361908e-01 1.537874e-01 1.728123e-01 1.933402e-01
#> [106] 2.153526e-01 2.387169e-01 2.632211e-01 2.886521e-01 3.148727e-01
#> [111] 3.418513e-01 3.696177e-01 3.981740e-01 4.274164e-01 4.571130e-01
#> [116] 4.869507e-01 5.166245e-01 5.459220e-01 5.747545e-01 6.031173e-01
#> [121] 6.310070e-01 6.583505e-01 6.849849e-01 7.106929e-01 7.352692e-01
#> [126] 7.585821e-01 7.805959e-01 8.013465e-01 8.208882e-01 8.392482e-01
#> [131] 8.564102e-01 8.723307e-01 8.869715e-01 9.003296e-01 9.124478e-01
#> [136] 9.234049e-01 9.332914e-01 9.421891e-01 9.501618e-01 9.572598e-01
#> [141] 9.635304e-01 9.690273e-01 9.738148e-01 9.779642e-01 9.815477e-01
#> [146] 9.846314e-01 9.872731e-01 9.895229e-01 9.914253e-01 9.930221e-01
#> [151] 9.943530e-01 9.954555e-01 9.963639e-01 9.971087e-01 9.977158e-01
#> [156] 9.982077e-01 9.986033e-01 9.989191e-01 9.991694e-01 9.993662e-01
#> [161] 9.995199e-01 9.996392e-01 9.997310e-01 9.998011e-01 9.998542e-01
#> [166] 9.998940e-01 9.999237e-01 9.999455e-01 9.999615e-01 9.999731e-01
#> [171] 9.999814e-01 9.999873e-01 9.999914e-01 9.999943e-01 9.999962e-01
#> [176] 9.999975e-01 9.999984e-01 9.999990e-01 9.999994e-01 9.999996e-01
#> [181] 9.999998e-01 9.999999e-01 9.999999e-01 1.000000e+00 1.000000e+00
#> [186] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [191] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [196] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [201] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [206] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [211] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [216] 1.000000e+00By design, similar to the ordinary DC-FFT algorithm by Biscarri, Zhao & Brunner (2018), its results are identical to the G-DC procedure, if \(n\) and the number of possible observed values is small. Thus, differences can be observed for larger numbers:
set.seed(1)
pp1 <- runif(250)
va1 <- sample(0:50, 250, TRUE)
vb1 <- sample(0:50, 250, TRUE)
pp2 <- pp1[1:248]
va2 <- va1[1:248]
vb2 <- vb1[1:248]
sum(abs(dgpbinom(NULL, pp1, va1, vb1, method = "DivideFFT")
- dgpbinom(NULL, pp1, va1, vb1, method = "Convolve")))
#> [1] 0
sum(abs(dgpbinom(NULL, pp2, va2, vb2, method = "DivideFFT")
- dgpbinom(NULL, pp2, va2, vb2, method = "Convolve")))
#> [1] 0The reason is that the G-DC-FFT method splits the input
probs, val_p and val_q vectors
into parts such that the numbers of possible observations of all parts
are as equally sized as possible. Their distributions are then computed
separately with the G-DC approach. The results of the portions are then
convoluted by means of the Fast Fourier Transformation. For small \(n\) and small distribution sizes, no
splitting is needed. In addition, the G-DC-FFT procedure, just like the
DC-FFT method, does not produce probabilities \(\leq 5.55e\text{-}17\), i.e. smaller values
are rounded off to \(0\), if the total
number of possible observations is smaller than \(750\), whereas the smallest possible result
of the DC algorithm is \(\sim
1e\text{-}323\). This is most likely caused by the used FFTW3
library.
The Generalized Discrete Fourier Transformation of the
Characteristic Function (G-DFT-CF) approach is requested with
method = "Characteristic".
set.seed(1)
pp <- runif(10)
wt <- sample(1:10, 10, TRUE)
va <- sample(0:10, 10, TRUE)
vb <- sample(0:10, 10, TRUE)
dgpbinom(NULL, pp, va, vb, wt, "Characteristic")
#> [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.837237e-16
#> [26] 6.270062e-16 1.364746e-15 2.935666e-15 6.201829e-15 1.292176e-14
#> [31] 2.657237e-14 5.394193e-14 1.081902e-13 2.144802e-13 4.201557e-13
#> [36] 8.135509e-13 1.557735e-12 2.949809e-12 5.527683e-12 1.025814e-11
#> [41] 1.885776e-11 3.434640e-11 6.196980e-11 1.106787e-10 1.956340e-10
#> [46] 3.425394e-10 5.948077e-10 1.025224e-09 1.753750e-09 2.972596e-09
#> [51] 4.985314e-09 8.275458e-09 1.362195e-08 2.227979e-08 3.622799e-08
#> [56] 5.845270e-08 9.332219e-08 1.473012e-07 2.302797e-07 3.576650e-07
#> [61] 5.529336e-07 8.496291e-07 1.292864e-06 1.943382e-06 2.888042e-06
#> [66] 4.257944e-06 6.248675e-06 9.128095e-06 1.322640e-05 1.893515e-05
#> [71] 2.675612e-05 3.741507e-05 5.199255e-05 7.194684e-05 9.895330e-05
#> [76] 1.347017e-04 1.809349e-04 2.399008e-04 3.150314e-04 4.112231e-04
#> [81] 5.341537e-04 6.888863e-04 8.788234e-04 1.106198e-03 1.374340e-03
#> [86] 1.690272e-03 2.065290e-03 2.511885e-03 3.037800e-03 3.641214e-03
#> [91] 4.311837e-03 5.039293e-03 5.824625e-03 6.686091e-03 7.651765e-03
#> [96] 8.740859e-03 9.945159e-03 1.122411e-02 1.252016e-02 1.378863e-02
#> [101] 1.502576e-02 1.627450e-02 1.759663e-02 1.902489e-02 2.052786e-02
#> [106] 2.201243e-02 2.336424e-02 2.450429e-02 2.543095e-02 2.622065e-02
#> [111] 2.697857e-02 2.776636e-02 2.855637e-02 2.924236e-02 2.969655e-02
#> [116] 2.983772e-02 2.967384e-02 2.929746e-02 2.883252e-02 2.836282e-02
#> [121] 2.788971e-02 2.734351e-02 2.663438e-02 2.570794e-02 2.457639e-02
#> [126] 2.331289e-02 2.201380e-02 2.075053e-02 1.954176e-02 1.836001e-02
#> [131] 1.716200e-02 1.592047e-02 1.464084e-02 1.335803e-02 1.211826e-02
#> [136] 1.095708e-02 9.886542e-03 8.897658e-03 7.972694e-03 7.098018e-03
#> [141] 6.270583e-03 5.496952e-03 4.787457e-03 4.149442e-03 3.583427e-03
#> [146] 3.083701e-03 2.641746e-03 2.249767e-03 1.902455e-03 1.596805e-03
#> [151] 1.330879e-03 1.102475e-03 9.084265e-04 7.447312e-04 6.071616e-04
#> [156] 4.918629e-04 3.956251e-04 3.158260e-04 2.502339e-04 1.968330e-04
#> [161] 1.537458e-04 1.192445e-04 9.179821e-05 7.010494e-05 5.308547e-05
#> [166] 3.984854e-05 2.965115e-05 2.187013e-05 1.598631e-05 1.157497e-05
#> [171] 8.295941e-06 5.881266e-06 4.121776e-06 2.854642e-06 1.953341e-06
#> [176] 1.320224e-06 8.809465e-07 5.799307e-07 3.763587e-07 2.406488e-07
#> [181] 1.515662e-07 9.401686e-08 5.742327e-08 3.451481e-08 2.039831e-08
#> [186] 1.184350e-08 6.751380e-09 3.777327e-09 2.073644e-09 1.116337e-09
#> [191] 5.887148e-10 3.036829e-10 1.529887e-10 7.516829e-11 3.598151e-11
#> [196] 1.676154e-11 7.585978e-12 3.326430e-12 1.407529e-12 5.717381e-13
#> [201] 2.216360e-13 8.149551e-14 2.825209e-14 9.182470e-15 2.781725e-15
#> [206] 7.813323e-16 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [211] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [216] 0.000000e+00
pgpbinom(NULL, pp, va, vb, wt, "Characteristic")
#> [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.837237e-16
#> [26] 9.107298e-16 2.275475e-15 5.211141e-15 1.141297e-14 2.433473e-14
#> [31] 5.090710e-14 1.048490e-13 2.130392e-13 4.275194e-13 8.476751e-13
#> [36] 1.661226e-12 3.218961e-12 6.168770e-12 1.169645e-11 2.195459e-11
#> [41] 4.081235e-11 7.515875e-11 1.371285e-10 2.478072e-10 4.434412e-10
#> [46] 7.859806e-10 1.380788e-09 2.406013e-09 4.159763e-09 7.132359e-09
#> [51] 1.211767e-08 2.039313e-08 3.401508e-08 5.629487e-08 9.252285e-08
#> [56] 1.509756e-07 2.442977e-07 3.915989e-07 6.218786e-07 9.795436e-07
#> [61] 1.532477e-06 2.382106e-06 3.674970e-06 5.618352e-06 8.506394e-06
#> [66] 1.276434e-05 1.901301e-05 2.814111e-05 4.136751e-05 6.030266e-05
#> [71] 8.705877e-05 1.244738e-04 1.764664e-04 2.484132e-04 3.473665e-04
#> [76] 4.820683e-04 6.630032e-04 9.029039e-04 1.217935e-03 1.629158e-03
#> [81] 2.163312e-03 2.852198e-03 3.731022e-03 4.837220e-03 6.211560e-03
#> [86] 7.901832e-03 9.967122e-03 1.247901e-02 1.551681e-02 1.915802e-02
#> [91] 2.346986e-02 2.850915e-02 3.433378e-02 4.101987e-02 4.867163e-02
#> [96] 5.741249e-02 6.735765e-02 7.858176e-02 9.110192e-02 1.048906e-01
#> [101] 1.199163e-01 1.361908e-01 1.537874e-01 1.728123e-01 1.933402e-01
#> [106] 2.153526e-01 2.387169e-01 2.632211e-01 2.886521e-01 3.148727e-01
#> [111] 3.418513e-01 3.696177e-01 3.981740e-01 4.274164e-01 4.571130e-01
#> [116] 4.869507e-01 5.166245e-01 5.459220e-01 5.747545e-01 6.031173e-01
#> [121] 6.310070e-01 6.583505e-01 6.849849e-01 7.106929e-01 7.352692e-01
#> [126] 7.585821e-01 7.805959e-01 8.013465e-01 8.208882e-01 8.392482e-01
#> [131] 8.564102e-01 8.723307e-01 8.869715e-01 9.003296e-01 9.124478e-01
#> [136] 9.234049e-01 9.332914e-01 9.421891e-01 9.501618e-01 9.572598e-01
#> [141] 9.635304e-01 9.690273e-01 9.738148e-01 9.779642e-01 9.815477e-01
#> [146] 9.846314e-01 9.872731e-01 9.895229e-01 9.914253e-01 9.930221e-01
#> [151] 9.943530e-01 9.954555e-01 9.963639e-01 9.971087e-01 9.977158e-01
#> [156] 9.982077e-01 9.986033e-01 9.989191e-01 9.991694e-01 9.993662e-01
#> [161] 9.995199e-01 9.996392e-01 9.997310e-01 9.998011e-01 9.998542e-01
#> [166] 9.998940e-01 9.999237e-01 9.999455e-01 9.999615e-01 9.999731e-01
#> [171] 9.999814e-01 9.999873e-01 9.999914e-01 9.999943e-01 9.999962e-01
#> [176] 9.999975e-01 9.999984e-01 9.999990e-01 9.999994e-01 9.999996e-01
#> [181] 9.999998e-01 9.999999e-01 9.999999e-01 1.000000e+00 1.000000e+00
#> [186] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [191] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [196] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [201] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [206] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [211] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [216] 1.000000e+00As can be seen, the G-DFT-CF procedure does not produce probabilities \(\leq 2.2e\text{-}16\), i.e. smaller values are rounded off to 0, most likely due to the used FFTW3 library.
To assess the performance of the exact procedures, we use the
microbenchmark package. Each algorithm has to calculate the
PMF repeatedly based on random probability and value vectors. The run
times are then summarized in a table that presents, among other
statistics, their minima, maxima and means. The following results were
recorded on an AMD Ryzen 9 5900X with 64 GiB of RAM and Manjaro KDE
Linux (26.0.2).
library(microbenchmark)
n <- 2500
set.seed(1)
va <- sample(1:50, n, TRUE)
vb <- sample(1:50, n, TRUE)
f1 <- function() dgpbinom(NULL, runif(n), va, vb, method = "DivideFFT")
f2 <- function() dgpbinom(NULL, runif(n), va, vb, method = "Convolve")
f3 <- function() dgpbinom(NULL, runif(n), va, vb, method = "Characteristic")
microbenchmark(f1(), f2(), f3(), times = 51)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> f1() 115.2568 115.9255 119.4541 116.1445 117.5055 244.6307 51
#> f2() 290.6293 291.8993 292.9224 292.7669 293.7977 295.9105 51
#> f3() 650.4590 710.5991 722.6389 726.8567 740.3368 753.4988 51Clearly, the G-DC-FFT procedure is the fastest one. It outperforms both the G-DC and G-DFT-CF approaches. The latter one needs a lot more time than the others. Generally, the computational speed advantage of the G-DC-FFT procedure increases with larger \(n\) (and \(m\)).