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+00
The 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+00
By design, as proposed by Biscarri, Zhao & Brunner (2018), its results are identical to the DC procedure, if n ≤ 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] 0
The 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 ≤ 750. In addition, the DC-FFT
procedure does not produce probabilities ≤ 5.55e-17, i.e. smaller values are
rounded off to 0, if n > 750, whereas the smallest
possible result of the DC algorithm is ∼ 1e-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+00
As can be seen, the DFT-CF procedure does not produce probabilities ≤ 2.22e-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+00
Obviously, the RF procedure does produce probabilities ≤ 5.55e-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 Windows 10 Education (22H2).
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() 55.33225 55.48742 55.85213 55.5676 55.74928 59.60287 51
#> f2() 126.77123 127.29649 127.86562 127.6447 128.04243 132.08101 51
#> f3() 167.96266 168.48935 169.02589 168.6778 169.31518 171.60450 51
#> f4() 91.62972 91.74005 92.68508 93.1054 93.20073 94.17543 51
Clearly, 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+00
The 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+00
By 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] 0
The 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 ≤ 5.55e-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 ∼ 1e-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+00
As can be seen, the G-DFT-CF procedure does not produce probabilities ≤ 2.2e-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 Windows 10
Education (22H2).
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() 326.4028 326.8275 331.2380 327.2105 328.7831 506.4497 51
#> f2() 776.7544 780.1921 782.0670 781.7714 784.0173 788.9935 51
#> f3() 613.5195 621.4151 624.7435 623.1188 625.0551 655.8785 51
Clearly, 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).