JimmyJim JimmyJim - 21 days ago 8
R Question

Ratio-of-Uniforms Distribution in R

I have an exercise, in which i have to create an algorithm as follows:

ratio of Uniforms is based on the fact that for a random variable X with density f(x) we can generate X from the desired density by calculating X = U/V for a pair (U, V ) uniformly distributed in the set

Af = {(u,v):0 < v ≤ f(u/v)}


Random points can be sampled uniformly in Af by rejection from the min- imal bounding rectangle, i.e., the smallest possible rectangle that contains Af .
It is given by (u−, u+) × (0, v+) where

v+ = max f(x), x

u− = minx f(x), x

u+ = maxx f(x)


Then the Ratio-of-Uniforms method consists of the following simple steps:


  1. Generate random number U uniformly in (u−, u+).

  2. Generate random number V uniformly in (0, v+).

  3. Set X ← U/V .

  4. If V 2 ≤ f(X) accept and return X.

  5. Else try again.



My code so far:

x <- cnorm(1, mean = 0, sd=1)

myrnorm <- function(pdf){
## call rou() n times
pdf <- function(x) {exp(-x^2/2)}
}
rou <- function(u, v) {
uplus <- 1
vplus <- 1
n <- 100
u <- runif(n, min=0, max=uplus)
v <- runif(n, min=0, max=vplus)
xi <- v/u
while(v < sqrt(xi)) {
if(v^2 <= xi)
return(xi)
}
}


myx <- myrnorm(1000)
hist(myx)


But I really dont know how to go on. Im ´lost with this exercise. I would be really grateful for any advise.

VMO VMO
Answer

Following example 1 in page 8 of this link and your sample code, I came up this solution:

ratioU <- function(nvals)
{
  h_x = function(x) exp(-x)
  # u- is b-, u+ is b+ and v+ is a in the example:
  uminus = 0
  uplus = 2/exp(1)
  vplus = 1
  X.vals <- NULL
  i <- 0
  repeat {
    i <- i+1
    u <- runif(1,0,vplus)
    v <- runif(1,uminus,uplus)
    X <- u/v
    if(v^2 <= h_x(X)) {
      tmp <- X
    }
    else {
      next
    }
    X.vals <- c(X.vals,tmp)
    if(length(X.vals) >= nvals) break
  }
  answer <- X.vals  
  answer
}

sol = ratioU(1000) 
par(mfrow=c(1,2))
hist(sol,breaks=50, main= "using ratioU",freq=F)
hist(rexp(1000),breaks = 50, main="using rexp from R",freq=F)
par(mfrow=c(1,1))

par(mfrow=c(1,2))
plot(density(sol))
plot(density(rexp(1000)))
par(mfrow=c(1,1))

A lot of the code may be optimized but I think it is good enough like this for this purpose. I hope this helps.

Comments