# Generating random samples from a custom distribution

I am trying to generate random samples from a custom pdf using R. My pdf is:

I generated uniform samples and then tried to transform it to my custom distribution. I did this by finding the cdf of my distribution ($F_{X}(x)$) and setting it to the uniform sample ($u$) and solving for $x$.

To generate a random sample with the above distribution, get a uniform sample $u \in[0,1]$ and solve for $x$ in

I implemented it in R and I don’t get the expected distribution. Can anyone point out the flaw in my understanding?

nsamples <- 1000;
x <- runif(nsamples);

f <- function(x, u) {
return(3/2*(x-x^3/3) - u);
}

z <- c();
for (i in 1:nsamples) {
# find the root within (0,1)
r <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root; z <- c(z, r); }  ## Answer It looks like you figured out that your code works, but @Aniko pointed out that you could improve its efficiency. Your biggest speed gain would probably come from pre-allocating memory for z so that you’re not growing it inside a loop. Something like z <- rep(NA, nsamples) should do the trick. You may get a small speed gain from using vapply() (which specifies the returned variable type) instead of an explicit loop (there’s a great SO question on the apply family). > nsamples <- 1E5 > x <- runif(nsamples) > f <- function(x, u) 1.5 * (x - (x^3) / 3) - u > z <- c() > > # original version > system.time({ + for (i in 1:nsamples) { + # find the root within (0,1) + r <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root
+   z <- c(z, r)
+ }
+ })
user  system elapsed
49.88    0.00   50.54
>
> # original version with pre-allocation
> z.pre <- rep(NA, nsamples)
> system.time({
+ for (i in 1:nsamples) {
+   # find the root within (0,1)
+   z.pre[i] <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root + } + }) user system elapsed 7.55 0.01 7.78 > > > > # my version with sapply > my.uniroot <- function(x) uniroot(f, c(0, 1), tol = 0.0001, u = x)$root
> system.time({
+   r <- vapply(x, my.uniroot, numeric(1))
+ })
user  system elapsed
6.61    0.02    6.74
>
> # same results

And you don’t need the ; at the end of each line (are you a MATLAB convert?).