Generating random numbers from a non-uniform distribution over a unit circle

1.2k Views Asked by At

I have a non-uniform distribution over a unit circle, peaking at $\frac{3}{\pi}$ at $(0,0)$ and tailing away to zero at the edge,

$f(x,y) = \frac{3}{\pi} \cdot (1 - \sqrt{x^2+y^2})$

I wish to generate a few random points within the unit circle according to this distribution.

Is there a straightforward method?

2

There are 2 best solutions below

0
On BEST ANSWER

Probably the easiest to implement would be rejection sampling. A simple implementation would be as follows.

Let $\mu = \max_{(x,y) \in [0,1]^2} f(x,y) = \frac3{\pi}$.

  1. Sample $x,\,y, \, u \sim \text{Unif}[0,1]$ independent uniform random variables.
  2. If $u > f(x,y)/\mu$ then return to step 1. Else accept $(x,y)$.

To see why this works you will want to read the wikipedia page cited above (or another source).


Example Implementation in R

library(dplyr)
library(ggplot2)

samples <- sapply(1:100000, function(i){
  success <- FALSE

  while(success == FALSE){
    x <- runif(1,-1,1)
    y <- runif(1,-1,1)
    u <- runif(1,0,1)

    if( u < 1- sqrt(x^2 + y^2) ){
      success <- TRUE
    }
  }

  return( c(x,y) )
}) %>% t %>% as.data.frame %>% select(x = V1, y = V2)

ggplot(samples, aes(x,y)) +  stat_density2d(aes(fill=..level..), geom="polygon")

enter image description here

0
On

An alternative solution, is to rephrase the distribution in Polar coordinates

$$r = \sqrt{x^2 + y^2}, \, \theta = \tan^{-1}\left(\frac{y}{x}\right).$$

In these coordinates the probability density function is given to be

$$\widetilde f(r, \theta) = \frac{3}{\pi}r(1-r),$$ where the additional factor of $r$ occurs because both $f,\widetilde f$ are defined by integrals and we have the change of variables $dx dy = r \, dr d\theta$.

This can be re-written as

$$ \widetilde f(r,\theta) = \left( \frac1{2\pi} \right) \times \bigg( 6 r\, (1-r) \bigg) = g(\theta) \, h(r)$$ where $g(\theta) = \frac1{2\pi}$ is the pdf of a uniform distribution on $[0,2\pi]$, and $h(r) = 6 r(1-r)$ is the pdf of a Beta(2,2) distribution on $[0,1]$.

In particular this means that $\theta, r$ are independent. So we can obtain samples from the distribution by independently sampling $$\theta \sim \text{Unif}[0,2 \pi], \qquad r \sim \text{Beta}(2,2).$$

Many software packages would include methods of sampling from a Beta distribution, but if this is not supported then as in the previous solution, rejection sampling could also be applied.