[BUG] Integration test shows that there are numerical problems for large dimension `d > 100`

[hsloot]

[BUG] Integration test shows that there are numerical problems for large dimension `d > 100`

[hsloot]

Some findings:

• The issue does not appear for low dimensions.
• The issue does seems to be present for all parametrisations for the intensities have to be approximated. In particular, it is not an issue for the constant, linear, or Poisson Bernstein function.
• Since the rescaled overall minimum is unit Exponential, this should not be an issue of the sample size.
• Some probabilities seem be lower that .Machine$double.eps (~1e-16) as seen in the reprex below. However, the difference between the sum over all intensities and the Bernstein function in d does increase but is not very high.
• There does not seem to be a similar problem for the Cuadras-Augé and the LFM samplers.

library(rmo)

ex_intensities_alpha_stable(d=125L, alpha=0.4)
#>   [1] 2.212888e-02 1.073568e-04 1.401013e-06 2.986698e-08 8.916098e-10
#>   [6] 3.434778e-11 1.616342e-12 9.051402e-14 5.908662e-15 4.399482e-16
#>  [11] 3.677403e-17 3.421603e-18 3.522201e-19 3.980770e-20 4.897729e-21
#>  [16] 6.515972e-22 9.337372e-23 1.437891e-23 2.373331e-24 4.183353e-25
#>  [21] 7.841928e-26 1.557942e-26 3.273315e-27 7.265646e-28 1.702352e-28
#>  [26] 4.204565e-29 1.092414e-29 2.978611e-30 8.505755e-31 2.540527e-31
#>  [31] 7.932439e-32 2.588649e-32 8.826091e-33 3.141466e-33 1.165745e-33
#>  [36] 4.503482e-34 1.808941e-34 7.549033e-35 3.272032e-35 1.472966e-35
#>  [41] 6.886524e-36 3.342832e-36 1.683727e-36 8.792549e-37 4.756490e-37
#>  [46] 2.663876e-37 1.544004e-37 9.260828e-38 5.748138e-38 3.692092e-38
#>  [51] 2.453664e-38 1.686616e-38 1.198671e-38 8.804383e-39 6.681780e-39
#>  [56] 5.238659e-39 4.242988e-39 3.550236e-39 3.068873e-39 2.740414e-39
#>  [61] 2.527704e-39 2.408015e-39 2.369053e-39 2.406888e-39 2.525286e-39
#>  [66] 2.736312e-39 3.062366e-39 3.540152e-39 4.227611e-39 3.316221e-39
#>  [71] 6.648331e-39 8.756859e-39 1.192016e-38 1.677216e-38 2.439799e-38
#>  [76] 3.670038e-38 5.710022e-38 9.191009e-38 1.530961e-37 2.639799e-37
#>  [81] 4.713339e-37 8.717613e-37 1.670893e-36 3.320148e-36 6.842264e-36
#>  [86] 1.463024e-35 3.247123e-35 7.484453e-35 1.792709e-34 4.465711e-34
#>  [91] 1.157970e-33 3.128449e-33 8.812942e-33 2.589982e-32 7.943072e-32
#>  [96] 2.542991e-31 8.505616e-31 2.976880e-30 1.092903e-29 4.221151e-29
#> [101] 1.719424e-28 7.394746e-28 3.355999e-27 1.605031e-26 8.085578e-26
#> [106] 4.300949e-25 2.431222e-24 1.474180e-23 9.671523e-23 6.888313e-22
#> [111] 5.306931e-21 4.393071e-20 3.898977e-19 3.747394e-18 3.997341e-17
#> [116] 4.873571e-16 6.882948e-15 1.113623e-13 2.040121e-12 4.353065e-11
#> [121] 1.177471e-09 4.365827e-08 2.306677e-06 2.197669e-04 3.467189e-01

^{Created on 2020-09-19 by the reprex package (v0.3.0)}

[hsloot]

Another finding: Normal qq-plots does not show a visual difference from theoretical normal quantiles ..

library(rmo)
library(ggplot2)

n <- 1e4L
d <- 125L
alpha <- 0.4
bf <- AlphaStableBernsteinFunction(alpha=alpha)

scale <- valueOf(bf, d, 0L)
x <- scale * apply(rmo_ex_arnold(n, d, ex_intensities_alpha_stable(d, alpha=alpha)), 1, min)

ks.test(x, pexp)
#> 
#>  One-sample Kolmogorov-Smirnov test
#> 
#> data:  x
#> D = 0.064275, p-value < 2.2e-16
#> alternative hypothesis: two-sided

ggplot(data.frame(x = qnorm(pexp(x))), aes(sample = x)) +
	geom_qq() + geom_qq_line()

^{Created on 2020-09-19 by the reprex package (v0.3.0)}

[hsloot]

In https://stat.ethz.ch/R-manual/R-devel/library/base/html/Random.html, it says

Do not rely on randomness of low-order bits from RNGs. Most of the supplied uniform generators return 32-bit integer
values that are converted to doubles, so they take at most 2^32 distinct values and long runs will return duplicated
values (Wichmann-Hill is the exception, and all give at least 30 varying bits.)

This explains at least why duplicate values occur suddenly in higher dimensions. This cannot be fixed and seems to be an aspect one has to accept, if working with R. An option could be to build interfaces for other 64bit RNG's.

[hsloot]

In the initializer of ExArnoldGenerator, the method binomial_coefficient is called. This function might induce undefined behaviour due to long integer overflow.

Possible solution: Create a recursion to fill the wt_generators_ and shock_generators_ vectors.