Statistical Testing of Random Number Generators in R, Exams of Latin language

Download Statistical Testing of Random Number Generators in R and more Exams Latin language in PDF only on Docsity! Package ‘randtoolbox’ April 19, 2009 Type Package Title toolbox for pseudo and quasi random number generation and RNG tests. Version 1.07 Date 2009-02-10 Author Christophe Dutang and Diethelm Wuertz (except some R port of the SFMT algorithm from M. Matsumoto, M. Saito and WELL generator from P. L’Ecuyer, Knuth-TAOCP RNG from D. Knuth and the Sobol sequence from the ACM transactions). See LICENCE file for details. Maintainer Christophe Dutang <christophe.dutang@ensimag.fr> Description The package provides (1) pseudo random generators - general linear congruential generators (Park Miller) and multiple recursive generators (Knuth TAOCP), generalized feedback shift register (SF-Mersenne Twister algorithm and WELL generator); (2) a quasi random generator - the Torus algorithm, the Sobol sequence, the Halton sequence (thus include Van der Corput sequence) and (3) some additional tests such as the gap test, the serial test, the poker test... For true random number generation, use the ’random’ package, for Latin Hypercube Sampling (a hybrid qmc method), use the ’lhs’ package, a number of RNGs and tests for RNGs are provided by ’RDieHarder’, all available on CRAN. Depends R (>= 2.6.0) License GPL (>= 2) Encoding latin1 LazyLoad yes Repository CRAN Date/Publication 2009-02-15 13:58:51 1 2 auxiliary R topics documented: auxiliary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 coll.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 freq.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 gap.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 order.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 poker.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 pseudo.randtoolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 quasi.randtoolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 serial.test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Index 23 auxiliary Auxiliary functions for ’randtoolbox’ package. Description Stirling numbers of the second kind and permutation of positive integers. Usage stirling(n) permut(n) Arguments n a positive integer. Details stirling computes stirling numbers of second kind i.e. Stirlkn = k ∗ Stirlkn−1 + Stirlk−1n−1 with Stirl1n = Stirl n n = 1. e.g. n = 0, returns 1 n = 1, returns a vector with 0,1 n = 2, returns a vector with 0,1,1 n = 3, returns a vector with 0,1,3,1 n = 4, returns a vector with 0,1,7,6,1... Go to wikipedia for more details. permut compute permutation of 1, ..., n and store it in a matrix. e.g. n = 1, returns a matrix with freq.test 5 References Planchet F., Jacquemin J. (2003), L’utilisation de methodes de simulation en assurance. Bulletin Francais d’Actuariat, vol. 6, 11, 3-69. (available online) L’Ecuyer P. (2001), Software for uniform random number generation distinguishing the good and the bad. Proceedings of the 2001 Winter Simulation Conference. (available online) L’Ecuyer P. (2007), Test U01: a C library for empirical testing of random number generators. ACM Trans. on Mathematical Software 33(4), 22. See Also other tests of this package freq.test, serial.test, poker.test, order.test and gap.test ks.test for the Kolmogorov Smirnov test and acf for the autocorrelation function. Examples # (1) poisson approximation # coll.test(runif) # (2) exact distribution # coll.test(SFMT, 2^7, 2^10, 10000) freq.test the Frequency test Description The Frequency test for testing random number generators. Usage freq.test(u, seq = 0:15, echo = TRUE) Arguments u sample of random numbers in ]0,1[. echo logical to plot detailed results, default TRUE seq a vector of contiguous integers, default 0:15. 6 freq.test Details We consider a vector u, realisation of i.i.d. uniform random variables U1, . . . , Un. The frequency test works on a serie seq of ordered contiguous integers (s1, . . . , sd), where sj ∈ ZZ. From the sample u, we compute observed integers as di = bui ∗ (sd + 1) + s1c, (i.e. di are uniformely distributed in {s1, . . . , sd}). The expected number of integers equals to j is m = 1sd−s1+1 × n. Finally, the chi-squared statistic is S = d∑ j=1 (card(di = sj)−m)2 m . Value a list with the following components : statistic the value of the chi-squared statistic. p.value the p-value of the test. observed the observed counts. expected the expected counts under the null hypothesis. residuals the Pearson residuals, (observed - expected) / sqrt(expected). Author(s) Christophe Dutang. References Planchet F., Jacquemin J. (2003), L’utilisation de methodes de simulation en assurance. Bulletin Francais d’Actuariat, vol. 6, 11, 3-69. (available online) L’Ecuyer P. (2001), Software for uniform random number generation distinguishing the good and the bad. Proceedings of the 2001 Winter Simulation Conference. (available online) L’Ecuyer P. (2007), Test U01: a C library for empirical testing of random number generators. ACM Trans. on Mathematical Software 33(4), 22. See Also other tests of this package gap.test, serial.test, poker.test, order.test and coll.test ks.test for the Kolmogorov Smirnov test and acf for the autocorrelation function. gap.test 7 Examples # (1) # freq.test(runif(1000)) print( freq.test( runif(1000000), echo=FALSE) ) # (2) # freq.test(runif(1000), 1:4) freq.test(runif(1000), 10:40) gap.test the Gap test Description The Gap test for testing random number generators. Usage gap.test(u, lower = 0, upper = 1/2, echo = TRUE) Arguments u sample of random numbers in ]0,1[. lower numeric for the lower bound, default 0. upper numeric for the upper bound, default 1/2. echo logical to plot detailed results, default TRUE Details We consider a vector u, realisation of i.i.d. uniform random variables U1, . . . , Un. The gap test works on the ’gap’ variables defined as Gi = { 1&if lower ≤ Ui ≤ upper\0&otherwise\ Let p the probability that Gi equals to one. Then we compute the length of zero gaps and denote by nj the number of zero gaps of length j. The chi-squared statistic is given by S = m∑ j=1 (nj − npj)2 npj , where pj stands for the probability the length of zero gaps equals to j ((1− p)2pj) and m the max number of lengths (at least ⌊ log(10−1)−2 log(1−p)−log(n) log(p) ⌋ ). 10 poker.test References Planchet F., Jacquemin J. (2003), L’utilisation de methodes de simulation en assurance. Bulletin Francais d’Actuariat, vol. 6, 11, 3-69. (available online) L’Ecuyer P. (2001), Software for uniform random number generation distinguishing the good and the bad. Proceedings of the 2001 Winter Simulation Conference. (available online) L’Ecuyer P. (2007), Test U01: a C library for empirical testing of random number generators. ACM Trans. on Mathematical Software 33(4), 22. See Also other tests of this package freq.test, serial.test, poker.test, gap.test and coll.test ks.test for the Kolmogorov Smirnov test and acf for the autocorrelation function. Examples # (1) mersenne twister vs torus # order.test(runif(6000)) order.test(torus(6000)) # (2) # order.test(runif(4000), 4) order.test(torus(4000), 4) # (3) # order.test(runif(5000), 5) order.test(torus(5000), 5) poker.test the Poker test Description The Poker test for testing random number generators. Usage poker.test(u , nbcard = 5, echo = TRUE) Arguments u sample of random numbers in ]0,1[. echo logical to plot detailed results, default TRUE nbcard a numeric for the number of cards, we assume that the length of u is a multiple of nbcard. poker.test 11 Details We consider a vector u, realisation of i.i.d. uniform random variables U1, . . . , Un. Let us note k the card number (i.e. nbcard). The poker test computes a serie of ’hands’ in {0, . . . , k − 1} from the sample hi = buidc (u must have a length dividable by k). Let nj be the number of ’hands’ with (exactly) j different cards. The probability is pj = k! kk(k − j)! ∗ Sjk ∗ ( j k )(k − j), where Sjk denotes the Stirling numbers of the second kind. Finally the chi-squared statistic is S = k−1∑ j=0 nj − npj/k)2 npj/k . Value a list with the following components : statistic the value of the chi-squared statistic. p.value the p-value of the test. observed the observed counts. expected the expected counts under the null hypothesis. residuals the Pearson residuals, (observed - expected) / sqrt(expected). Author(s) Christophe Dutang. References Planchet F., Jacquemin J. (2003), L’utilisation de methodes de simulation en assurance. Bulletin Francais d’Actuariat, vol. 6, 11, 3-69. (available online) L’Ecuyer P. (2001), Software for uniform random number generation distinguishing the good and the bad. Proceedings of the 2001 Winter Simulation Conference. (available online) L’Ecuyer P. (2007), Test U01: a C library for empirical testing of random number generators. ACM Trans. on Mathematical Software 33(4), 22. See Also other tests of this package freq.test, serial.test, gap.test, order.test and coll.test ks.test for the Kolmogorov Smirnov test and acf for the autocorrelation function. 12 pseudo.randtoolbox Examples # (1) hands of 5 'cards' # poker.test(runif(50000)) # (2) hands of 4 'cards' # poker.test(runif(40000), 4) # (3) hands of 42 'cards' # poker.test(runif(420000), 42) pseudo.randtoolbox Toolbox for pseudo and quasi random number generation Description General linear congruential generators such as Park Miller sequence, generalized feedback shift register such as SF-Mersenne Twister algorithm and WELL generator; and a quasi random generator (pseudo random generators) and the Torus algorithm (quasi random generation). Usage congruRand(n, dim = 1, mod = 2^31-1, mult = 16807, incr = 0, echo) SFMT(n, dim = 1, mexp = 19937, usepset = TRUE, withtorus = FALSE, usetime = FALSE) WELL(n, dim = 1, order = 512, temper = FALSE, version = "a") knuthTAOCP(n, dim = 1) setSeed(seed) Arguments n number of observations. If length(n) > 1, the length is taken to be the required number. dim dimension of observations (must be <=100 000, default 1). seed a single value, interpreted as a positive integer for the seed. e.g. append your day, your month and your year of birth. mod an integer defining the modulus of the linear congruential generator. mult an integer defining the multiplier of the linear congruential generator. incr an integer defining the increment of the linear congruential generator. echo a logical to plot the seed while computing the sequence. mexp an integer for the mersenne exponent of SFMT algorithm. see details withtorus a numeric in ]0,1] defining the proportion of the torus sequence appended to the SFMT sequence; or a logical equals to FALSE (default). pseudo.randtoolbox 15 See Also .Random.seed for what is done in R about random number generation. Examples # (1) the Park Miller sequence # # Park Miller sequence, i.e. mod = 2^31-1, mult = 16807, incr=0 # the first 10 seeds used in Park Miller sequence # 16807 1 # 282475249 2 # 1622650073 3 # 984943658 4 # 1144108930 5 # 470211272 6 # 101027544 7 # 1457850878 8 # 1458777923 9 # 2007237709 10 setSeed(1) congruRand(10, echo=TRUE) # the 9998+ th terms # 925166085 9998 # 1484786315 9999 # 1043618065 10000 # 1589873406 10001 # 2010798668 10002 setSeed(1614852353) #seed for the 9997th term congruRand(5, echo=TRUE) # (2) the SF Mersenne Twister algorithm SFMT(1000) #Kolmogorov Smirnov test #KS statistic should be around 0.037 ks.test(SFMT(1000), punif) #KS statistic should be around 0.0076 ks.test(SFMT(10000), punif) #different mersenne exponent with a fixed parameter set # SFMT(10, mexp = 607, usepset = FALSE) SFMT(10, mexp = 1279, usepset = FALSE) SFMT(10, mexp = 2281, usepset = FALSE) SFMT(10, mexp = 4253, usepset = FALSE) SFMT(10, mexp = 11213, usepset = FALSE) SFMT(10, mexp = 19937, usepset = FALSE) SFMT(10, mexp = 44497, usepset = FALSE) SFMT(10, mexp = 86243, usepset = FALSE) 16 pseudo.randtoolbox SFMT(10, mexp = 132049, usepset = FALSE) SFMT(10, mexp = 216091, usepset = FALSE) #use different sets of parameters [default when possible] # for(i in 1:7) print(SFMT(1, mexp = 607)) for(i in 1:7) print(SFMT(1, mexp = 2281)) for(i in 1:7) print(SFMT(1, mexp = 4253)) for(i in 1:7) print(SFMT(1, mexp = 11213)) for(i in 1:7) print(SFMT(1, mexp = 19937)) #use a fixed set and a fixed seed #should be the same output setSeed(08082008) SFMT(1, usepset = FALSE) setSeed(08082008) SFMT(1, usepset = FALSE) # (3) withtorus argument # # one third of outputs comes from Torus algorithm u <- SFMT(1000, with=1/3) # the third term of the following code is the first term of torus sequence print(u[666:670] ) # (4) WELL generator # # 'basic' calls # WELL512 WELL(10, order = 512) # WELL1024 WELL(10, order = 1024) # WELL19937 WELL(10, order = 19937) # WELL44497 WELL(10, order = 44497) # WELL19937 with tempering WELL(10, order = 19937, temper = TRUE) # WELL44497 with tempering WELL(10, order = 44497, temper = TRUE) # tempering vs no tempering setSeed(08082008) WELL(10, order =19937) setSeed(08082008) WELL(10, order =19937, temper=TRUE) # (5) Knuth TAOCP generator # knuthTAOCP(10) knuthTAOCP(10, 2) pseudo.randtoolbox 17 # (6) How to set the seed? # all example is duplicated to ensure setSeed works # congruRand setSeed(1302) congruRand(1) setSeed(1302) congruRand(1) # SFMT setSeed(1302) SFMT(1, usepset=FALSE) setSeed(1302) SFMT(1, usepset=FALSE) # BEWARE if you do not set usepset to FALSE setSeed(1302) SFMT(1) setSeed(1302) SFMT(1) # WELL setSeed(1302) WELL(1) setSeed(1302) WELL(1) # Knuth TAOCP setSeed(1302) knuthTAOCP(1) setSeed(1302) knuthTAOCP(1) # (7) computation times on my macbook, mean of 1000 runs # ## Not run: # algorithm time in seconds for n=10^6 # classical Mersenne Twister 0.066 # SF Mersenne Twister 0.044 # WELL generator 0.065 # Knuth TAOCP 0.046 # Park Miller 0.108 n <- 1e+06 mean( replicate( 1000, system.time( runif(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( SFMT(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( WELL(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( knuthTAOCP(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( congruRand(n), gcFirst=TRUE)[3]) ) ## End(Not run) 20 quasi.randtoolbox setSeed(6) torus(5) #the same setSeed(1) torus(10) #no use of the machine time torus(10, use=FALSE) #Kolmogorov Smirnov test #KS statistic should be around 0.0019 ks.test(torus(1000), punif) #KS statistic should be around 0.0003 ks.test(torus(10000), punif) #the mixed Torus sequence torus(10, mix=TRUE) par(mfrow = c(1,2)) acf(torus(10^6)) acf(torus(10^6, mix=TRUE)) # (2) Halton sequences # # uniform variate halton(n = 10, dim = 5) # normal variate halton(n = 10, dim = 5, normal = TRUE) # some plots par(mfrow = c(2, 2), cex = 0.75) hist(halton(n = 5000, dim = 1), main = "Uniform Halton", xlab = "x", col = "steelblue3", border = "white") hist(halton(n = 5000, dim = 1), main = "Normal Halton", xlab = "x", col = "steelblue3", border = "white") # (3) Sobol sequences # # uniform variate sobol(n = 10, dim = 5, scrambling = 3) # normal variate sobol(n = 10, dim = 5, scrambling = 3, normal = TRUE) # some plots hist(sobol(5000, 1, scrambling = 2), main = "Uniform Sobol", xlab = "x", col = "steelblue3", border = "white") hist(sobol(5000, 1, scrambling = 2, normal = TRUE), main = "Normal Sobol", serial.test 21 xlab = "x", col = "steelblue3", border = "white") # (4) computation times on my macbook, mean of 1000 runs # ## Not run: # algorithm time in seconds for n=10^6 # Torus algo 0.058 # mixed Torus algo 0.087 # Halton sequence 0.878 # Sobol sequence 0.214 n <- 1e+06 mean( replicate( 1000, system.time( torus(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( torus(n, mixed=TRUE), gcFirst=T)[3]) ) mean( replicate( 1000, system.time( halton(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( sobol(n), gcFirst=TRUE)[3]) ) ## End(Not run) serial.test the Serial test Description The Serial test for testing random number generators. Usage serial.test(u , d = 8, echo = TRUE) Arguments u sample of random numbers in ]0,1[. echo logical to plot detailed results, default TRUE d a numeric for the dimension, see details. When necessary we assume that d is a multiple of the length of u. Details We consider a vector u, realisation of i.i.d. uniform random variables U1, . . . , Un. The serial test computes a serie of integer pairs (pi, pi+1) from the sample u with pi = buidc (u must have an even length). Let nj be the number of pairs such that j = pi × d + pi+1. If d=2, we count the number of pairs equals to 00, 01, 10 and 11. Since all the combination of two elements in {0, . . . , d− 1} are equiprobable, the chi-squared statistic is S = d−1∑ j=0 nj − n/(2d2))2 n/(2d2) . 22 serial.test Value a list with the following components : statistic the value of the chi-squared statistic. p.value the p-value of the test. observed the observed counts. expected the expected counts under the null hypothesis. residuals the Pearson residuals, (observed - expected) / sqrt(expected). Author(s) Christophe Dutang. References Planchet F., Jacquemin J. (2003), L’utilisation de methodes de simulation en assurance. Bulletin Francais d’Actuariat, vol. 6, 11, 3-69. (available online) L’Ecuyer P. (2001), Software for uniform random number generation distinguishing the good and the bad. Proceedings of the 2001 Winter Simulation Conference. (available online) L’Ecuyer P. (2007), Test U01: a C library for empirical testing of random number generators. ACM Trans. on Mathematical Software 33(4), 22. See Also other tests of this package freq.test, gap.test, poker.test, order.test and coll.test ks.test for the Kolmogorov Smirnov test and acf for the autocorrelation function. Examples # (1) # serial.test(runif(1000)) print( serial.test( runif(1000000), d=2, e=FALSE) ) # (2) # serial.test(runif(5000), 5)