Rellene aleatoriamente una cuadrícula 3D según una función de densidad de probabilidad p (x, y, z)

¿Cómo puedo llenar una cuadrícula 3D en el orden especificado por una función de densidad de probabilidad dada?

Usando python, me gustaría establecer puntos en un orden aleatorio , pero de acuerdo con alguna distribución de probabilidad especificada sobre esa región, sin puntos repetidos.

Secuencialmente:

  • crear una cuadrícula 3D discreta
  • especifique una función de densidad de probabilidad para cada punto de la cuadrícula, pdf (x, y, z)
  • establecer un punto (x0, y0, z0) cuya ubicación aleatoria es proporcional al pdf (x, y, z)
  • Continúe agregando puntos (sin repetir) hasta que se llenen todas las ubicaciones

El resultado deseado es una lista de todos los puntos (sin repeticiones) de todos los puntos de la cuadrícula, para que se rellenen.

Lo siguiente no implementa el dibujo de un gaussiano multivariado:

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

Esto se debe a que p(x)*p(y)*p(z) != p(x,y,z) menos que las tres variables sean independientes. Puede considerar algo como un muestreador de Gibbs para extraer de la distribución conjunta mediante un dibujo secuencial de distribuciones univariadas.

En el caso específico de la multivariable normal, puede usar (ejemplo completo)

 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) 

Aquí hay un ejemplo, usando un pdf gaussiano (ver gráfico). Este código se adapta fácilmente a cualquier pdf especificado:

 %matplotlib qt import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #number of points to lay down: n = 4000; #create meshgrid: min, max, L = -5, 5, 91; [x_grid,y_grid,z_grid] = np.meshgrid(np.linspace(min,max,L),np.linspace(min,max,L),np.linspace(min,max,L)) xi,yi,zi = x_grid.ravel(),y_grid.ravel(),z_grid.ravel() #create normalized pdf (gaussian here): pdf = np.exp(-(x_grid**2 + y_grid**2 + z_grid**2)); pdf = pdf/np.sum(pdf); #obtain indices of randomly selected points, as specified by pdf: randices = np.random.choice(np.arange(x_grid.ravel().shape[0]), n, replace = False,p = pdf.ravel()) #random positions: x_rand = xi[randices] y_rand = yi[randices] z_rand = zi[randices] fig = plt.figure(); ax = fig.add_subplot(111, projection='3d',aspect='equal') svals = 16; ax.scatter(x_rand, y_rand, z_rand, s=svals, alpha=.1) 

gráfico de dispersión generado por código