Ratio-of-Uniforms Distribution in R -
i have exercise, in have create algorithm follows:
ratio of uniforms based on fact random variable x density f(x) can generate x desired density calculating x = u/v pair (u, v ) uniformly distributed in set
af = {(u,v):0 < v ≤ f(u/v)}
random points can sampled uniformly in af rejection min- imal bounding rectangle, i.e., smallest possible rectangle contains af . given (u−, u+) × (0, v+) where
v+ = max f(x), x u− = minx f(x), x u+ = maxx f(x)
then ratio-of-uniforms method consists of 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 , return x.
else try again.
my code 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 dont know how go on. im ´lost exercise. grateful advise.
following example 1 in page 8 of link , sample code, came solution:
ratiou <- function(nvals) { h_x = function(x) exp(-x) # u- b-, u+ b+ , v+ in example: uminus = 0 uplus = 2/exp(1) vplus = 1 x.vals <- null <- 0 repeat { <- 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 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 code may optimized think enough this purpose. hope helps.
Comments
Post a Comment