19
$\begingroup$

I am trying to generate random samples from a custom pdf using R. My pdf is: $$f_{X}(x) = \frac{3}{2} (1-x^2), 0 \le x \le 1$$

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$.

$$ F_{X}(x) = \Pr[X \le x] = \int_{0}^{x} \frac{3}{2} (1-y^2) dy = \frac{3}{2} (x - \frac{x^3}{3}) $$

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);
}
$\endgroup$
7
  • 1
    $\begingroup$ Must be a coding mistake. I don't use R, so I can't say what the mistake is exactly -- but I just coded up your solution (taking care to take the middle root of the cubic polynomial, which always lies between 0 and 1), and I get good agreement between the samples and the expected distribution. Could it be a problem with your root finder? What's wrong with the samples you're getting? $\endgroup$
    – jpillow
    Commented Jul 9, 2011 at 21:58
  • $\begingroup$ I tried your code (which is not very efficient, by the way) and do get the expected distribution. $\endgroup$
    – Aniko
    Commented Jul 9, 2011 at 22:20
  • $\begingroup$ @jpillow and @Aniko My mistake. When I used nsamples <- 1e6 it was a good match. $\endgroup$
    – Anand
    Commented Jul 9, 2011 at 22:32
  • 2
    $\begingroup$ @Anand One way is to observe that $x = 2\sin(\arcsin(u)/3)$, allowing direct calculation of $x$ in terms of $u$. $\endgroup$
    – whuber
    Commented Jul 11, 2011 at 14:08
  • 1
    $\begingroup$ @Anand en.wikipedia.org/wiki/… $\endgroup$
    – whuber
    Commented Jul 14, 2011 at 16:08

1 Answer 1

11
$\begingroup$

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?).

$\endgroup$
2
  • $\begingroup$ Thanks for your detailed answer and for pointing out vapply. I have been coding in C/C++ for a very long time and that is the reason for ; affliction! $\endgroup$
    – Anand
    Commented Jul 10, 2011 at 3:16
  • 1
    $\begingroup$ +1 It looks like replacing uniroot by a direct solution speeds things up by three more orders of magnitude: you should have no trouble obtaining around $10^7$ variates per second. $\endgroup$
    – whuber
    Commented Jul 11, 2011 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.