how to generate a set of random points within a half-circle uniformly in python or R?

Question

I am trying to generate a set of points within a half-circle uniformly.

size = 1000
t  = np.random.random(size)*np.pi*2
u  = np.random.random(size) + np.random.random(size)
r = np.where(u > 1, 2 - u, u)
x = r*cos(t)
y = r*sin(t)
coor = (x,y)
for idx, value in enumerate(y):
    if value<0:
        x[idx]=-3
        y[idx]=-3
f, ax = plt.subplots(figsize = (3,3))
plt.scatter(x, y)

this piece of code has 2 bugs.

1. the plot has a lot (-3,-3) points to be removed
2. it seems to be inefficient.

the figure shown as below is NOT uniformly, as the points are more of center than others.

another plot shown as below could be viewed as uniformly.

any idea to fix the bugs would be appreciated.

Answer 1

You should generate uniformly distributed angles phi, and take the sqrt of uniformly generated radius r (which takes into account that we want to sample uniformly on the area, see explanation below), to make sure you sample points uniformly in the half-circle.

import numpy as np
import matplotlib.pyplot as plt

# sample
size = 10000
R = 1
phi = np.random.random(size=size) * np.pi
r = np.sqrt(np.random.random(size=size)) * R

# transform
x = r * np.cos(phi)
y = r * np.sin(phi)

# plot
f = plt.figure(figsize=(12,12))
a = f.add_subplot(111)
a.scatter(x, y, marker='.')
a.set_aspect('equal')
plt.show()

Explanation

To generate uniformly distributed points on a (half-)circle we have to make sure that each infinitesimal area or segment is "hit" with the same probability. We can simply sample phi from a uniform random distribution [0, 1), multiplied by np.pi (so [0, pi)), since all angles should have the same probability to be sampled. But if we sample r from a uniform random distribution in [0, 1), we generate too many points at small radia, and not enough at large radia, since the area grows like r**2. To take that fact into account, we have to bias our sampled radia accordingly, and in this case we can apply the biasing by simply taking the square root (np.sqrt), to apply the correct weighting to the sampled radius values, and take into account the larger area of the outer rings.

A much better and more thorough explanation is found here: https://stackoverflow.com/a/50746409/1170207

Performance comparison to rejection sampling methods

Since this method is basically an inversion sampling method, we compare its performance to a rejection sampling algorithm.

import numpy as np
x, y = np.random.random(size=(2,10000))
%timeit r, phi = np.sqrt(x), y
# 19 µs ± 33.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit m = x**2 + y**2 <= 1; xx, yy = x[m], y[m]
# 81.5 µs ± 271 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

With the rejection sampling method we also cannot ensure that we draw a select number of variates, so we have to repeat the process until we have. This cannot be vectorized that nicely, unless we accept to sample too many values, and discard additional ones.

Answer 2

You should generate points in the enclosing rectangle, and remove points which are not in the half-circle

# generate about n points in a half-circle with
# given center (x, y) and given radius, and y>0

points <- function(x, y, radius, n) {
    n2 = n * 4 / pi # each point has pi/4 probability to survive

    # make [-1, 1] * [-1, 1] square
    xs = runif(n2, -1, 1)
    ys = runif(n2, 0, 1)  # or just runif(n2)
    points = cbind(xs, ys)

    # keep only points in circle with center (0,0) and radius 1 with y>0
    ind = (xs**2 + ys**2 <= 1) # the condition ys>=0 is obeyed already
    points = points[ind,]

    # move/stretch to given center and radius
    points = points * radius
    points[,1] = points[,1] + x
    points[,2] = points[,2] + y
}

# generate about 1000 points with center(1,1) and radius 3
points = f(1, 1, 3, 1000)

# plot them, making them smaller for better visibility
plot(points, cex=0.3)

Answer 3

In R:

runif_in_semicircle <- function(n, radius=1){
  theta <- runif(n, 0, pi)
  r <- radius * sqrt(runif(n))
  cbind(r*cos(theta), r*sin(theta))
}

sims <- runif_in_semicircle(1000)
plot(sims[,1], sims[,2], asp=1, pch=19)

We can check it works by evaluating an integral.

# integrand example
f <- function(x) x[1]^2 + exp(x[2])

set.seed(666)
sims <- runif_in_semicircle(10000)
fsims <- apply(sims, 1, f)
mean(fsims)*pi/2 # approximates the integral of f over the half-disk
# 2.890905

Now we numerically evaluate the integral of f.

library(SphericalCubature)
adaptIntegrateBallPolar(f, n=2, lowerLimit = 0, upperLimit = pi)
# $integral
# [1] 2.880598

Answer 4

You can try out following in R:

size <- 1000
maxRad <- 1 #maximum radius of the half-circle
r <- runif(1000,0,maxRad) #generate random radius
phi <- runif(size,0,pi) #generate angle for polarcoordinates (between 0 and pi since you want a halfcircle)
x <- r*cos(phi) #polarcoordinates
y <- r*sin(phi)
plot(x,y)

You can put it in a function

halfCircle <- function(size, maxRad) {
 r <- runif(1000,0,maxRad)
 phi <- runif(size,0,1)*pi
 x <- r*cos(phi)
 y <- r*sin(phi)
 plot(x,y)
}

and try out whether it delivers you "more acceptable random" results.

Answer 5

This method uses Python, and a rejection algorithm.

First create uniformly distributed variables in a square. Then reject all points that lie outside the half circle you are interested in. Here I chose the simplest case, but it can be extended to any half circle if you add rotation, scaling and translation of these points.

The code is as follows:

import numpy as np
import matplotlib.pyplot as plt

n = 10000
r = 3  # radius

uniform_square = np.random.uniform(-r,r,(int(2.65*n),2))
radius = np.sqrt(uniform_square[:,0]**2 + uniform_square[:,1]**2)

uniform_circle = uniform_square[radius<=r,:]

uniform_half_circle = uniform_circle[uniform_circle[:,0]>=0,:]

final_points = uniform_half_circle[0:n,:]

fig, axs = plt.subplots(1, 1)
axs.scatter(final_points[:,0], final_points[:,1], marker='.')
axs.set_aspect('equal', 'box')
plt.show()

Since it uses numpy for all steps, its relatively quick. To get a fixed number of points, generate more than you need initially and downsize the array. Since I started with a perfect square (area=4) and only need a half circle (area = pi/2) we need to generate about 2.6 times as many points as we need in the end.

The for loop in python is not very fast, so try to stick to numpy only functions and operations.

Answer 6

Not sure what you mean 'uniformly'.

Here is one approach to generate 'uniformly' distributed points along x- and y-axis, but not very efficient

import numpy as np
import matplotlib.pyplot as plt

x = np.random.random(size) * 2 - 1
y = np.random.random(size)
r = np.sqrt(x**2 + y**2)
x[r > 1] = -3
y[r > 1] = -3
plt.plot(x, y, 'o')
plt.show()