JimmyJim - 21 days ago 8

R Question

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:

- Generate random number U uniformly in (u−, u+).
- Generate random number V uniformly in (0, v+).
- Set X ← U/V .
- If V 2 ≤ f(X) accept and return X.
- 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.

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.

Source (Stackoverflow)

Comments