[R] uniform generator (default)

[Alec Stephenson]

Recieving digests.
> RNGkind(NULL)
[1] "Marsaglia-Multicarry" "Kinderman-Ramage"

I would appreciate it if anybody has any comments on the following.
Please do not comment on the R functions themselves, since they merely
mimic a (bivariate simplification of a) C routine called from S.

In particular, I would like to know if anything is available with regard
to the theoretical properties of the default uniform generator under
transformations such as these. 

testing <-
function(n,alpha) {
  sim <- matrix(0,ncol=2,nrow=n)
  for(i in 1:n) 
  { 
    u <- runif(1)
    # runif(1)  ####
    if(runif(1) < alpha) z <- -log(runif(1))-log(runif(1))
    else z <- -log(runif(1))
    sim[i,1] <- exp(-(z * u^alpha))
    sim[i,2] <- exp(-(z * (1-u)^alpha))
  }
  sim
}

tmp <- testing(100000,0.5) # about 30secs 4 me
plot(tmp[,1],tmp[,2],pch=".")

testing <-
function(n,alpha) {
  sim <- matrix(0,ncol=2,nrow=n)
  for(i in 1:n) 
  { 
    u <- runif(1)
    runif(1)  ####
    if(runif(1) < alpha) z <- -log(runif(1))-log(runif(1))
    else z <- -log(runif(1))
    sim[i,1] <- exp(-(z * u^alpha))
    sim[i,2] <- exp(-(z * (1-u)^alpha))
  }
  sim
}

tmp <- testing(100000,0.5) # about 30secs 4 me
plot(tmp[,1],tmp[,2],pch=".")


> version
         _                
platform i686-pc-linux-gnu
arch     i686             
os       linux-gnu        
system   i686, linux-gnu  
status                    
major    1                
minor    3.1              
year     2001             
month    08               
day      31               
language R  

As an aside, the routine written specifically in R may take the form

"testing.r"<- 
function(n = 1, alpha = 1, margins = "frechet")
{
        if(length(n) != 1 || mode(n) != "numeric" || n <= 0 || n%%1 != 0)
            stop("`n' must be a positive integer")
        if(length(alpha) != 1 || mode(alpha) != "numeric" || alpha <= 0 ||
        alpha > 1) stop("`alpha' must be a parameter in (0,1]")
        ####
	u <- runif(n)
        z <- rgamma(n, 1+rbinom(n,1,alpha))
	sim <- z * cbind(x=u, y=1-u)^alpha
        exp(-sim)
} 

which doesn't have the same structured output (for obvious reasons).

Incedently, the distribution function from which these sample from is the
copula F(x,y)=exp(-(u^(1/alpha)+v^(1/alpha))^alpha) where u=-log(x) and
v=-log(y).

Thanks in advance.

Alec Stephenson                            tel +44 (0) 1524 593950
Department of Mathematics and Statistics   fax +44 (0) 1524 592681
Lancaster University            email a.stephenson at lancaster.ac.uk
Lancaster LA1 4YF



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._