[R] How to generate the random numbers uniformly distributed on the unit disc?

[Scionforbai]

It depends on what you're actually doing, but I normally prefer Van
der Corput pseudo-random sequences to draw uniformly AND regularly
distributed points on a 2d disc or 3d sphere.
VdC sequences need binary and ternary decompositions of the numbers of
simulated points. They are a little bit time consuming, but the
sequences can be

Here a small piece of code (all functions of mine begin with M.):

#######################################################

# binary or ternary decomposition. results in a vector.
M.decomp <- function(n,base) {
   l <- trunc(log(n,base)+1)
   dec <- vector(mode = "integer", length = l)

   for (i in c(1:l) ) {
      dec[l-i+1] <- n%%base
      n <- n%/%base
   }
   return(dec)

}

# Calculates the i-th terms of VdC binary and ternary sequences:
M.seq.vdc <- function(i) {
   # binary decomp.
   v2 <- M.decomp(i,2)
   lv <- v2/(2**(length(v2):1))
   # ternary decomp.
   u3 <- M.decomp(i,3)
   lu <- u3/(3**(length(u3):1))

   # actual terms of VdC sequence:
   vdc2=sum(lv)
   vdc3=sum(lu)

   return(c(vdc2,vdc3))
}

# Uniformly distributed points on the unitary disc (2D)
M.vdc.2d <- function(n,plot=F) {

   a <- lapply(1:n,M.seq.vdc)
   a <- matrix(unlist(a),ncol=2,byrow=T)
   u <- 2*pi*a[,1]
   v <- sqrt(a[,2])
   x <- v*cos(u)
   y <- v*sin(u)

   if (plot) plot(x,y,pch=".",asp=1)

   invisible(list(phi = u, rho = v, x = x, y = y))
}

# For comparison, the "uniform" algorithm
M.unif.disc <- function (n = 10^4, R = 1,plot=F) {
   phi <- runif(n) * 2 * pi
   rho <- R * sqrt(runif(n))
   xx <- rho * cos(phi)
   yy <- rho * sin(phi)

   if (plot) plot(xx, yy, pch = ".", asp = 1)

   invisible(list(phi = phi, rho = rho, x = xx, y = yy))
}


# Uniformly distributed points on the unitary sphere (3D)
M.vdc.3d <- function(n) {

   a <- lapply(1:n,M.seq.vdc)
   out   <- matrix(unlist(a),ncol=2,byrow=T)
   u <- out[,1]
   v <- out[,2]
   x <- cos(2*pi*u)*sqrt(1-v^2)
   y <- sin(2*pi*u)*sqrt(1-v^2)
   out <- cbind(x,y,v)
   return(out)
}


####################################################"


# test and visualization:
op <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
par("mar" = c(1, 1, 1, 1) + 0.1)

vdc <- M.vdc.2d(n,plot=T)
uni <- M.unif.disc(n,plot=T)
par(op)

Actually, one should randomly rotate the VdC points (and maybe
permutate the realization?) for each simulation (the sequence is
deterministic). Google for "Van der Corput sequence" to find more
details about its nice properties.

Regards,

Marco