It is possible to generate points within an ellipse without using rejection sampling too by carefully considering its definition in polar form. From wikipedia the polar form of an ellipse is given by
Intuitively speaking, we should sample polar angle θ more often where the radius is larger. Put more mathematically, our PDF for the random variable θ should be p(θ) dθ = dA / A, where dA is the area of a single segment at angle θ with width dθ. Using the equation for polar angle area dA = 1/2 r2 dθ and the area of an ellipse being π a b, then the PDF becomes
To randomly sample from this PDF, one direct method is the inverse CDF technique. This requires calculating the cumulative density function (CDF) and then inverting this function. Using Wolfram Alpha to get the indefinite integral, then inverting it gives inverse CDF of
where u runs between 0 and 1. So to sample a random angle θ, you just generate a uniform random number u between 0 and 1, and substitute it into this equation for the inverse CDF.
To get the random radius, the same technique that works for a circle can be used (see for example Generate a random point within a circle (uniformly)).
Here is some sample Python code which implements this algorithm:
import numpy
import matplotlib.pyplot as plt
import random
# Returns theta in [-pi/2, 3pi/2]
def generate_theta(a, b):
u = random.random() / 4.0
theta = numpy.arctan(b/a * numpy.tan(2*numpy.pi*u))
v = random.random()
if v < 0.25:
return theta
elif v < 0.5:
return numpy.pi - theta
elif v < 0.75:
return numpy.pi + theta
else:
return -theta
def radius(a, b, theta):
return a * b / numpy.sqrt((b*numpy.cos(theta))**2 + (a*numpy.sin(theta))**2)
def random_point(a, b):
random_theta = generate_theta(a, b)
max_radius = radius(a, b, random_theta)
random_radius = max_radius * numpy.sqrt(random.random())
return numpy.array([
random_radius * numpy.cos(random_theta),
random_radius * numpy.sin(random_theta)
])
a = 2
b = 1
points = numpy.array([random_point(a, b) for _ in range(2000)])
plt.scatter(points[:,0], points[:,1])
plt.show()