I draw the attached plot drawn with the following code. Each color represents a cell. How can I find the average number of cells that one user (any point within the square) can be in? By average I mean "the user can be anywhere within the whole area". For example, averaging over 10000 random locations of a point within the square area.

enter image description here

 ell = t [Function] t^-20.75];500, 2}]; 
 H = RandomVariate[ExponentialDistribution[1], 50]; 
 r = Table[H[[i]] ell[Norm[{x, y} - X[[i]]]], {i, 1, 20}]; 
 s = Table[r[[i]]/( Total[Delete[r, i]] + 30.99), {i, 1, 20}]; 

 RegionPlot[s[[i]] >= 0.1, {x, -0.7, 0.7}, {y, -0.7, 0.7}, 
 MaxRecursion -> 2, PlotPoints -> 40, 
 PlotStyle -> {FaceForm[{Opacity[0.5], ColorData[97][i]}]}, 
 BoundaryStyle -> ColorData[97][i]], {i, 1, 20}], Graphics[Point[X]]]
  • $\begingroup$ Wouldn't it be easier (and more correct) to sum up the areas of the shaded regions (with multiplicity) and to divide the result by the area of the square? $\endgroup$ Jul 30 '20 at 10:27
  • $\begingroup$ @HenrikSchumacher, seems to be a very good idea. But I do not know how to do it. I really appreciate if you can help me on this. $\endgroup$ Jul 30 '20 at 10:39
  • 1
    $\begingroup$ Sum[Area@ImplicitRegion[f >= 0.1, {{x, -0.7, 0.7}, {y, -0.7, 0.7}}],{f, s}]/(2 0.7)^2? See also the documentation of Area for a detailed explanation of the options of Area that allow you to control the accuracy of the computation (AccuracyGoal and PrecisionGoal). $\endgroup$ Jul 30 '20 at 10:45
  • $\begingroup$ @HenrikSchumacher, Thanks a lot. It seems to be computationally very expensive approach. Is there any alternative way that you can think of!! $\endgroup$ Jul 30 '20 at 11:14
  • 1
    $\begingroup$ Why do you think so? It is not more expensive to compute an approximate area of an impicit region than computing an approximate triangle mesh via RegionPlot. Basically that is what Area does: It meshes the domain and sums the triangle areas. $\endgroup$ Jul 30 '20 at 11:17

Well, a Monte-Carlo approach could look like this. There is a considerable amount of number crunching involved, so we better compile the working horse function:

cf = Compile[{{H, _Real, 1}, {X1, _Real, 1}, {X2, _Real, 1}, {Z, _Real, 1}},
   Block[{n, R, r, x1, x2},
    x1 = Compile`GetElement[Z, 1];
    x2 = Compile`GetElement[Z, 2];
    r = H ((X1 - x1)^2 + (X2 - x2)^2)^-1.875;
    R = Total[r] + 3.99;
    Total[UnitStep[r - (0.1/1.1) R]]
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"


pts = RandomReal[{-0.7, 0.7}, {10000000, 2}];
{X1, X2} = Transpose[X];
N@Mean@cf[H, X1, X2, pts]


The function cf counts how many entries of the vector s are above 0.1 for a given {x,y} = Z. Applying cf in a listable way to all elements of pts and taking the mean yields the empirical average number of satisfied inequalities.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.