JimmyJim - 1 year ago 82
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.

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
}
}

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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download