rou.R

X <- numeric(length=1000)
params.a <- sqrt(30) / 4
params.c <- 0.811

sqrt.fx <- function(x) {
  f <- (x ^ 4) - 2 * (x ^ 3) + (x ^ 2)
  return(sqrt(30 * f))
}

sample <- function() {
  while(TRUE) {
    u <- runif(1, min=0, max=params.a)
    v <- runif(1, min=0, max=params.c)
    
    ratio <- v / u
    if (u <= sqrt.fx(ratio)) {
      break
    }
  }
  return(ratio)
}

for (i in 1:1000) {
  X[i] <- sample()
}
hist(X)
# plot(density(X))