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

fX(x)=32(1−x2),0≤x≤1I generated uniform samples and then tried to transform it to my custom distribution. I did this by finding the cdf of my distribution (FX(x)) and setting it to the uniform sample (u) and solving for x.

FX(x)=Pr

To generate a random sample with the above distribution, get a uniform sample u \in[0,1] and solve for x in \frac{3}{2} (x – \frac{x^3}{3}) = u

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
> head(z)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738
> head(z.pre)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738
> head(r)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738
```

And you don’t need the `;`

at the end of each line (are you a MATLAB convert?).

**Attribution***Source : Link , Question Author : Anand , Answer Author : Community*