mardi 28 février 2023

Generate a Uniform Spherical Distribution using rejection methods in Python

I've been trying to generate a uniform spherical distribution in Python using uniform random sampling. For some reason my presumed spherical distribution looks more like an ovoid than like a sphere. I am using the fact that a sphere is defined: x^2 + y^2 + z^2 = R^2 and assuming R = 1 I get that the condition that the points should satisfy to be inside the sphere is: x^2 + y^2 + z^2 <= 1. For some reason this does not work. I have a perfect circular distribution from a top view (projection in the xy plane), but there is clearly a elliptical geometry in the planes xz and yz. ``` import numpy as np import numpy.random as rand import matplotlib.pyplot as plt

N = 10000

def sample_cube(number_of_trials):
    """
    This function receives an integer and returns a uniform sample of points in a square.
    The return object is a list of lists, which consists of three entries, each is a list
    of the copordinates of a point in the space.
    """

    cube = rand.uniform(low = -1, high = 1, size = (number_of_trials, 3))
    x_cube = cube[:,0]
    y_cube = cube[:,1]
    z_cube = cube[:,2]
    return [x_cube, y_cube, z_cube]

def sample_sphere(cube):
    """
    This function takes a list on the form [x_cube, y_cube, z_cube] and then sample
    an spherical distribution of points.
    """
    in_sphere = np.where(np.sqrt(cube[0]**2 + cube[1]**2 + cube[2]**2) <= 1)
    return in_sphere

"""
Main Code
"""
cube = sample_cube(N)
sphere = sample_sphere(cube)
print(sphere[0])
print(cube)

print(cube[0][sphere[0]])
x_in_sphere = cube[0][sphere[0]]
y_in_sphere = cube[1][sphere[0]]
z_in_sphere = cube[2][sphere[0]]
fig = plt.figure()
ax = plt.axes(projection = "3d")
ax.scatter(x_in_sphere, y_in_sphere, z_in_sphere, s = 1, color = "black")
plt.show()
plt.clf
![Distribution generated by the code above](https://i.stack.imgur.com/N83Ui.png)
![Side view of the distribution (plane xz)](https://i.stack.imgur.com/Di5kJ.png)
![Top view of distribution (plane xy)](https://i.stack.imgur.com/8uQvD.png)

I was simply trying to get a uniform sphere. There should be something wrong with the approach, but I can not spot the mistake.



Aucun commentaire:

Enregistrer un commentaire