Test occupancy in arbitrary region - Monte Carlo Integration - Uniformity test

Hello all,

I want to do a Monte Carlo integration (because integral limits make it nontrivial) and of course the sampling method is important in such instances, as visually demonstrated by John D'Errico in a related question:

https://nl.mathworks.com/matlabcentral/answers/442977-monte-carlo-integration-in-matlab-help

My own question on the matter is: given an arbitrary region (maybe even 3d or higher, where the visualisation is hard or impossible), is there a metric or a test to verify if the sampling is reasonably uniform over the space sampled? For example if someone takes the code from the related question, is it possible to obtain a global measure of the points distribution uniformity?

An image to think about it: 2 x 10^5 points generated inside an arbitrary alphaShape, I want to check if points are indeed uniformly distributed, but visualisation does not help much here. Searching for more I found this is related to a Complete Spatial Randomness test if I understood correctly.

[William Rose]

@Gustavo Lunardon,

If I were doing this:

I would try to find a way to make an N-dimensional array (where N is the dimesnion of the region of integration). Call it pop, for population. The array elements ill be (hyper)cubes of equal volume, and together they will cover the region.

The initial value of the array elements that are outside the region or which straddle the border is is -1. THe initial value of all cubes that are totally inside the region will be zero. This is a weak point of my idea, because it may be difficult to efficiently determine which cubes are totally inside the region. (One approach: check each of the 2^N corner points of each cube.. If they are all inside then the cube is totally inside. That could perhaps occasionally fail for concave regions, but probably unlikely.)

As you do the integration by choosing random points, add 1 to the appropriate element of pop, each time you choose a point. When you are done, you reshape pop into a 1D vector and you throw away all elements whose vaue is -1. This values in the cells of this vector should be uniformly distributed, if your MOnte Calro sampled the space randomly and evenly. You can use the Kolmogorov-Smirnov test or the chi-squared test to test if the values are consistent with the uniform distribution.

Good luck!

[Gustavo Lunardon]

Hello @William Rose,

In fact, I had already tried something similar. This would allow an analytical approach to the problem but it did not work (left it running for a lot of time and no answer).

I had based my code on this MATLAB post, which is somewhat similar to what you proposed:

https://nl.mathworks.com/matlabcentral/answers/510273-integral-evaluation-in-an-alphashape

However, I have a matrix of functions I need to compute and within 25min it did not compute even the first term under the 'iterated' method so I tried changing my approach to a Monte Carlo integration. So now I am looking for some way of judging if the sampling method is OK, be it a single number, a plot, or something of the sort.

[Torsten]

Does the region you want to integrate over have the same dimension as the underlying space ? Or is it a hyperplane or of even lower dimension ?

[Gustavo Lunardon]

Hello @Torsten. The alphaShape in my case limits the nonconvex region of the integral, so in the 3d case it is the volume I am integrating over. I am not particularly interested in a plane that crosses it.

[Torsten]

The usual way to get a uniform distribution is to choose a hypercube of volume V that encloses your alphaShape and generate N random points uniformly distributed within the hypercube. Then select those n points that are within the alphaShape and compute the volume of the alphaShape V' as approximately

V' ~ V*n/N

or for the integral I over a function f as

I ~ V'/n * sum_{i=1}^{i=n} f(x_i).

[Gustavo Lunardon]

Thank you for your comment. I am following the rejection sampling method like you suggested. I think it should be OK then, if I generate points uniforming distributed and 'cut' a region from that space the points inside the region should be uniformly distributed as well, I suppose.

For the volume of the alphaShape there is a command volume(shp), where shp is the alphaShape object.

Thanks for confirming the formula

Is there a difference between your small letter n and capital N?

[Torsten]

n instead of N ; you only use the points within the alphaShape:

I ~ V'/n * sum_{i=1}^{i=n} f(x_i)

I'd first do tests for regions and functions you know the answer of in order to get a feeling for how many sampling points you need to achieve a sufficient accuracy.

[Gustavo Lunardon]

I see. In my counter I only increase N when it is valid so in my code N = n (I generated points inside a while instead of inside a for loop). Thank you for the clarification! Indeed, yesterday I generated points inside an alphaShape that was a rectangle and integrated for the simple limits symbollically and using the Monte Carlo approach (sample size 2 x 10^5) and got virtually the same result, so I think the code is working well now.

[Torsten]

You will also have to keep track of N together with n in order to compute V'.

Or to insert it in

I ~ V/N * sum_{i=1}^{i=n} f(x_i).

[Gustavo Lunardon]

Normally yes but for the specific case of alphaShape I can get the volume inside the polyhedron by the command

  % https://nl.mathworks.com/help/matlab/ref/alphashape.volume.html
  V_prime = volume(shp)

Instead of doing V' ~ V*n/N (which works even in case there is no alphaShape available). Do you see any problem in doing that?

[Torsten]

Do you see any problem in doing that?

No, of course not.

[Gustavo Lunardon]

Great. Thank you for your help! Happy holidays :)