La forma más fácil de lograr lo que está buscando es probablemente la función de integración de scipy. Aquí tu ejemplo:

 from scipy import integrate # Make data. def func(x,y,z): return x**y**z ranges = [[0,1], [0,1], [0,1]] result, error = integrate.nquad(func, ranges)

¿Eres consciente de que la función que creaste es diferente a la que muestras en la imagen? El que creaste es un exponencial (x^y^z) mientras que el que estás mostrando son solo multiplicaciones. Si quieres representar la función en la imagen, usa

 def func(x,y,z): return x*y*z

Espero que esto responda a tu pregunta, de lo contrario, ¡solo escribe un comentario!

Editar:

Leí mal tu publicación. Si solo tiene los resultados y no están espaciados regularmente, tendría que encontrar alguna forma de interpolación (es decir, lineal) y una tabla de búsqueda. Si no sabes cómo crear eso, házmelo saber. El resto de la respuesta indicada aún podría usarse si define func para devolver valores interpolados de sus datos originales

Una buena manera de hacerlo sería usar la integración scipy de tplquad . Sin embargo, para usar eso, necesitamos una función y no un punto de nube.

Una forma fácil de evitarlo es usar un interpolador, para obtener una función que se aproxime a nuestro punto de nube; por ejemplo, podemos usar RegularGridInterpolator de RegularGridInterpolator si los datos están en una cuadrícula regular:

 import numpy as np from scipy import integrate from scipy.interpolate import RegularGridInterpolator # Make data. def function(x,y,z): return x*y*z N = 5 xmin, xmax = 0, 1 ymin, ymax = 0, 1 zmin, zmax = 0, 1 x = np.linspace(xmin, xmax, N) y = np.linspace(ymin, ymax, N) z = np.linspace(zmin, zmax, N) values = function(*np.meshgrid(x,y,z, indexing='ij')) # Interpolate: function_interpolated = RegularGridInterpolator((x, y, z), values) # tplquad integrates func(z,y,x) f = lambda z,y,x : my_interpolating_function([z,y,x]) result, error = integrate.tplquad(f, xmin, xmax, lambda _: ymin, lambda _:ymax,lambda *_: zmin, lambda *_: zmax)

En el ejemplo anterior, obtenemos result = 0.12499999999999999 - ¡lo suficientemente cerca!

Según sus requisitos, parece que la técnica más adecuada sería la integración de Monte Carlo:

 # Step 0 start with some empirical data observed_points = np.random.uniform(0,1,size=(10000,3)) unknown_fn = lambda x: np.prod(x) # just used to generate fake values observed_values = np.apply_along_axis(unknown_fn, 1, observed_points) K = 1000000 # Step 1 - assume that f(x,y,z) can be approximated by an interpolation # of the data we have (you could get really fancy with the # selection of interpolation method - we'll stick with straight lines here) from scipy.interpolate import LinearNDInterpolator f_interpolate = LinearNDInterpolator(observed_points, observed_values) # Step 2 randomly sample from within convex hull of observed data # Step 2a - Uniformly sample from bounding 3D-box of data lower_bounds = observed_points.min(axis=0) upper_bounds = observed_points.max(axis=0) sampled_points = np.random.uniform(lower_bounds, upper_bounds,size=(K, 3)) # Step 2b - Reject points outside of convex hull... # Luckily, we get a np.nan from LinearNDInterpolator in this case sampled_values = f_interpolate(sampled_points) rejected_idxs = np.argwhere(np.isnan(sampled_values)) # Step 2c - Remember accepted values of estimated f(x_i, y_i, z_i) final_sampled_values = np.delete(sampled_values, rejected_idxs, axis=0) # Step 3 - Calculate estimate of volume of observed data domain # Since we sampled uniformly from the convex hull of data domain, # each point was selected with P(x,y,z)= 1 / Volume of convex hull volume = scipy.spatial.ConvexHull(observed_points).volume # Step 4 - Multiply estimated volume of domain by average sampled value I_hat = volume * final_sampled_values.mean() print(I_hat)

Para una derivación de por qué esto funciona, vea esto: https://cs.dartmouth.edu/wjarosz/publications/dissertation/appendixA.pdf

La primera respuesta explica muy bien el enfoque principal para manejar esto. Solo quería ilustrar una forma alternativa mostrando el poder del paquete sklearn y la regresión del aprendizaje automático.

Hacer la malla en 3D da una matriz numpy muy grande,

 import numpy as np N = 5 xmin, xmax = 0, 1 ymin, ymax = 0, 1 zmin, zmax = 0, 1 x = np.linspace(xmin, xmax, N) y = np.linspace(ymin, ymax, N) z = np.linspace(zmin, zmax, N) grid = np.array(np.meshgrid(x,y,z, indexing='ij')) grid.shape = (3, 5, 5, 5) # 2*5*5*5 = 250 numbers

Lo cual visualmente no es muy intuitivo con 250 números. Con diferente indexación posible ('ij' o 'xy'). Usando la regresión podemos obtener el mismo resultado con pocos puntos de entrada (15-20).

 # building random combinations from (x,y,z) X = np.random.choice(x, 20)[:,None] Y = np.random.choice(y, 20)[:,None] Z = np.random.choice(z, 20)[:,None] xyz = np.concatenate((X,Y,Z), axis = 1) data = np.multiply.reduce(xyz, axis = 1)

Entonces, la entrada (cuadrícula) es solo una matriz numérica 2D,

 xyz.shape (20, 3)

Con los datos correspondientes,

 data.shape = (20,)

Ahora la función de regresión e integración,

 from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LinearRegression from sklearn.pipeline import Pipeline from scipy import integrate pipe=Pipeline([('polynomial',PolynomialFeatures(degree=3)),('modal',LinearRegression())]) pipe.fit(xyz, data) def func(x,y,z): return pipe.predict([[x, y, z]]) ranges = [[0,1], [0,1], [0,1]] result, error = integrate.nquad(func, ranges)

 print(result) 0.1257

Este enfoque es útil con un número limitado de puntos.