Randomly fill a 3D grid according to a probability density function p(x,y,z)
Question: **How can I fill a 3D grid in the order specified by a given probability
density function?**

Using python, I'd like to lay down points in a _random_ order, but according
to some specified probability distribution over that region, with no repeated
points.

Sequentially:

  * create a discrete 3D grid
  * specify a probability density function for every grid point, pdf(x,y,z)
  * lay down a point (x0,y0,z0) whose random location is proportional to the pdf(x,y,z)
  * continue adding additional points, not recording entries if a location has already been filled
  * repeat until all spaces are filled

The desired result is a list of all points (no repeats) of all the points in
the grid, in order that they were filled.


Answer: The below does not implement drawing from a multivariate gaussian:

    
    
    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())
    

That is because `p(x)*p(y)*p(z) != p(x,y,z)` unless the three variables are
independent. You can consider something like a [Gibbs
sampler](https://en.wikipedia.org/wiki/Gibbs_sampling) to draw from the joint
distribution by sequentially drawing from univariate distributions.

In the specific case of the multivariate normal, you can use (full example)

    
    
    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 *
    
    num_points = 4000
    
    sigma = .5;
    mean = [0, 0, 0]
    cov = [[sigma**2,0,0],[0,sigma**2,0],[0,0,sigma**2]]
    
    x,y,z = np.random.multivariate_normal(mean,cov,num_points).T
    
    svals = 16
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d',aspect='equal')
    ax.scatter(x,y,z, s=svals, alpha=.1,cmap=cm.gray)
    

