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.74048 19.81483 20.15701 19.87764 19.97383 25.92402 51
#> f2() 49.08046 49.39233 49.62233 49.50199 49.75656 52.02728 51
#> f3() 123.69639 124.03941 124.72059 124.24852 124.50770 135.37129 51
#> f4() 140.28894 143.56082 144.45182 144.52264 145.17602 147.60842 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() 110.0922 110.5917 113.6268 110.8240 111.5613 240.0725 51
#> f2() 290.3807 292.0037 292.7153 292.5672 293.1404 297.4572 51
#> f3() 645.7101 707.3148 718.9250 722.5912 736.3879 756.5021 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\)).