Laying points on a grid, ordered by a probability density function (pdf)

Using python, I'd like to fill in/sort grid of points in random order. Specifically: how do I order/sort points in 3D according to a specified probability density function (pdf)?

The basic logic of what I'm trying to do in steps:

• a probability density function is defined for every grid point, p(x,y,z)
• fill a first position (x0,y0,z0) according to the pdf p(x,y,z)
• continue to fill additional points (xi,yi,zi) according to p(x,y,z), ignoring those points that overlap with previously filled points
• repeat until all spaces are filled
• result is an ordering of unique points pi = (xi,yi,zi) in order that they were filled

Code:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
from math import *

%matplotlib qt 

#saturated
num_points = 12000
#begins to show square shape
num_points = 4000
#looks ok?
num_points = 500

L = 3;
dl = .2;
sigma = .5;
[x_grid,y_grid,z_grid] = np.meshgrid(np.arange(-L,L,dl),np.arange(-L,L,dl),np.arange(-L,L,dl))
pdf = np.exp(-(x_grid**2 + y_grid**2 + z_grid**2)/sigma**2)/np.sum(np.exp(-(x_grid**2 + y_grid**2 + z_grid**2)/sigma**2))

xi_sorted = np.random.choice(x_grid.ravel(),x_grid.ravel().shape, replace=False, p = pdf.ravel())
yi_sorted = np.random.choice(x_grid.ravel(),x_grid.ravel().shape, replace=False, p = pdf.ravel())
zi_sorted = np.random.choice(x_grid.ravel(),x_grid.ravel().shape, replace=False, p = pdf.ravel())

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d',aspect='equal')
#ax.set_aspect('equal')

svals = 16;

ax.scatter(xi_sorted[0:num_points], yi_sorted[0:num_points], zi_sorted[0:num_points], s=svals, alpha=.1,cmap=cm.gray)

For a gaussian pdf, these points should fill in spherically symmetrically, but my method shows a distinctly cubic shape well before they hit the boundary of x = 3. Thanks!