# Ratio of uniform random numbers

# 1. Normal distribution

a1 <- 0   # Setting the enclosing rectangle
b1 <- 1
a2 <- -sqrt(2*exp(-1))
b2 <- -a2
# computing the region G for the pair (u,v) in polar coordinates
b <- tan(seq(-1.57,1.57,0.001)) # slopes for a regular grid of angles
u <- exp(-0.25*b^2) # u-coordinate of the boundary of G on the line v=b*u
plot(u,b*u,type="l",xlim=c(a1,b1),ylim=c(a2,b2),xaxs="i",yaxs="i",xlab="u",
     ylab="v",main="ratio of uniforms for the normal distribution")
abline(h=0)
# Illustration: The probability 0.5 < X < 0.55 where X=V/U
abline(0,0.5,col=2) #event V/U < 0.5
abline(0,0.55,col=2) # event V/U < 0.55
abline(v=exp(-1/16),col=2) # vertical line at u = exp(-0.25 * (0.5)^2) 
0.0005*sum(u*u*(1+b*b))/((b1-a1)*(b2-a2))
     #approximate acceptance probability (Riemann sum in polar coordinates)

f.myrnorm <- function(nrep)
{
  ## Purpose: normal random numbers by the ratio of uniforms
  ## -------------------------------------------------------------------------
  ## Arguments: nrep=number of replicates
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 17 Dec 2001, 14:22
  x <- rep(0,nrep)
  n <- 0
  a <- sqrt(2*exp(-1))
  while (n < nrep) {
    u <- runif(1)
    v <- a*(1-2*runif(1))
    if (v^2 < - 4*u^2*log(u)) {
      n <- n+1
      x[n] <- v/u
    }
  }
  x
}

qqnorm(f.myrnorm(1000))

#2. Gamma(g,1)-distribution (same procedure as in the normal case)
g <- 10.5
a1 <- 0
b1 <- exp(0.5*(g-1)*(log(g-1)-1))
a2 <- 0
b2 <- exp(0.5*(g+1)*(log(g+1)-1))
b <- tan(seq(0.001,1.57,0.001)) 
u <- exp(0.5*((g-1)*log(b) - b)) 
plot(u,b*u,type="l",xlim=c(a1,b1),ylim=c(a2,b2),xaxs="i",yaxs="i",xlab="u",
     ylab="v",main="ratio of uniforms for gamma distribution")
0.0005*sum(u*u*(1+b*b))/(b1*b2)