QUnif: Quasi Randum Numbers via Halton Sequences

In sfsmisc: Utilities from 'Seminar fuer Statistik' ETH Zurich

Quasi Randum Numbers via Halton Sequences

Description

These functions provide quasi random numbers or space filling or low discrepancy sequences in the p-dimensional unit cube.

Usage

sHalton(n.max, n.min = 1, base = 2, leap = 1)
 QUnif (n, min = 0, max = 1, n.min = 1, p, leap = 1, silent = FALSE)

Arguments

n.max	maximal (sequence) number.
n.min	minimal sequence number.
n	number of p-dimensional points generated in QUnif. By default, n.min = 1, leap = 1 and the maximal sequence number is n.max = n.min + (n-1)*leap.
base	integer \ge 2: The base with respect to which the Halton sequence is built.
min, max	lower and upper limits of the univariate intervals. Must be of length 1 or p.
p	dimensionality of space (the unit cube) in which points are generated.
leap	integer indicating (if > 1) if the series should be leaped, i.e., only every leapth entry should be taken.
silent	logical asking to suppress the message about enlarging the prime table for large p.

Value

sHalton(n,m) returns a numeric vector of length n-m+1 of values in [0,1].

QUnif(n, min, max, n.min, p=p) generates n-n.min+1 p-dimensional points in [min,max]^p returning a numeric matrix with p columns.

Note

For leap Kocis and Whiten recommend values of L=31,61,149,409, and particularly the L=409 for dimensions up to 400.

Author(s)

Martin Maechler

References

James Gentle (1998) Random Number Generation and Monte Carlo Simulation; sec.\ 6.3. Springer.

Kocis, L. and Whiten, W.J. (1997) Computational Investigations of Low-Discrepancy Sequences. ACM Transactions of Mathematical Software 23, 2, 266–294.

Examples

32*sHalton(20, base=2)

stopifnot(sHalton(20, base=3, leap=2) ==
          sHalton(20, base=3)[1+2*(0:9)])
## ------- a 2D Visualization -------

Uplot <- function(xy, axes=FALSE, xlab="", ylab="", ...) {
  plot(xy, xaxs="i", yaxs="i", xlim=0:1, ylim=0:1, xpd = FALSE,
       axes=axes, xlab=xlab, ylab=ylab, ...)
  box(lty=2, col="gray40")
}

do4 <- function(n, ...) {
  op <- mult.fig(4, main=paste("n =", n,": Quasi vs. (Pseudo) Random"),
                 marP=c(-2,-2,-1,0))$old.par
  on.exit(par(op))
  for(i in 1:2) {
     Uplot(QUnif(n, p=2), main="QUnif", ...)
     Uplot(cbind(runif(n), runif(n)), main="runif", ...)
  }
}
do4(100)
do4(500)
do4(1000, cex = 0.8, col="slateblue")
do4(10000, pch= ".", col="slateblue")
do4(40000, pch= ".", col="slateblue")

sfsmisc documentation built on Sept. 11, 2024, 6:53 p.m.