How to perform a double (numerical) integration of $f(x,y)$ over an irregular sample of $x$ and $y$ values I have arrays of $x$ data and $y$ data which can be thought of as being randomly distributed (the distribution is not actually random but the relationship between x and y is unknown). 
I also have a function $f(x.y)$ which I wish to integrate over $x_{min} \rightarrow x_{max}$ and $y_{min} \rightarrow y_{max}$.
I know that for a single integral of, say $f(x)$, a sample of $x$ data, the trapezium rule generalises nicely to irregular values of x as the the width of the trapeziums = $\Delta x = x_{i+1} - x_i $.  
The trapezium rule also generalises to regularly spaced data for double integrals, but what about double intgrals over irregularly spaced data samples?
Any suggestions would be most welocome,  I'm sure there must be a standard method to this problem but have not been able to find one! 
(PS. I am coding this in python fyi)
 A: I know that Voronoi diagrams are used for interpolating random points data. First you need to create the partitioning, for example using Fortune's algorithm.
The rest depends on the accuracy you need. If you only want an estimation, you can use piecewise constant interpolation (i.e. the value for each cell is the same as the value for the point you built it around). 
In this case the integral is just the sum of volumes $\sum_k f(x_k,y_k) S_k$ where $S_k$ is the area of the Voronoi cell around the point $x_k,y_k$.
If you want a more precise measurement, you need to interpolate properly. One good method is Natural_neighbor_interpolation.

There's also a dirty trick I like for estimation.
Let's say you have an array of $f_k=f(x_k,y_k)$ defined at randomly distributed points.
Then use the Gaussian sum for interpolation with some $a$ of the order of the average distance between the points:
$$f(x,y) \approx \sum_k A_k \exp \left(-\frac{(x-x_k)^2}{a^2}-\frac{(y-y_k)^2}{a^2}\right)$$
The coefficients are found from a system of linear equations:
$$\sum_k A_k \exp \left(-\frac{(x_n-x_k)^2}{a^2}-\frac{(y_n-y_k)^2}{a^2}\right)=f_n$$
Because Gaussian function falls fast, the integral then can be approximated as the integral over the whole space:
$$I \approx \pi a^2 \sum_k A_k$$
Of course, this method is not accurate, but if you choose a good value for $a$, you can at least do as good as piecewise constant interpolation, and without the need to compute Voronoi diagrams.
The above assumes the distances between the points are all of the same order, i.e. the points are evenly distributed. Otherwise, you would need to use several values of $a$ to account for sparser and denser areas.
