Problem with extracting specific contour from a ContourPlot

I am trying to exctract contours of specific value (2.6*10^-6) out of a ContourPlot and measure their area.

Below you can see the code. It is noticable that not all of the relevant contours showed in the left plot are shown on the right one as well. Any suggestion of what am I doing wrong?

data = Import["C:\\Users\\doa\\Documents\\private\\school\\PhD\\thesistex\\sigma_limit_difference.xlsx", {"Sheets", "synthetic clusters"}];
data = data[[2 ;; 3002, 2 ;; 3]];
\[ScriptCapitalD] = SmoothKernelDistribution[data, 9.34];
contour = ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0,
1024}, PlotRange -> {Automatic}, PlotPoints -> 50, Contours -> 10, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.7]];
contour2 =ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0,
1024}, PlotRange -> {Automatic}, PlotPoints -> 60, Contours -> {2.6*10^-6}, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.7],ContourLabels -> True];
{contour, contour2}


I am not sure if there is a need to insert the excel file (in case it is needed please explain how and I will)

thanks in advance

Doron

All thank you very much for the help:

I tried again, this time using some random data so everyone can folow up.

The solution of changing the PlotRange parameter to "All" did'nt solve the problem. It is not just the "White areas" are seen on the left plot with plotrange full, but most impoprtant is whay dont we see them all in the right Plot? Maybe it is related to very small contours?

As for the contour area: Now I get an answer (Thanks...) but still get this warning message:

"NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>"

here is the code I am using:

SeedRandom[5]; data = RandomReal[1024, {3000, 2}];
\[ScriptCapitalD] = SmoothKernelDistribution[data, 9.34];
contour = ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0,
1024}, PlotRange -> {All}, PlotPoints -> 50, Contours -> 10, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.4]];
contour2 =ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y,0,1024}, PlotRange -> {All}, PlotPoints -> 60, Contours -> {2.52*10^-6}, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.4], ContourLabels -> True];
{contour, contour2}
GG = NIntegrate[Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}]

-
The white indicates that the function's value is not in the (automatically selected) plot range. Try PlotRange -> All and see if it looks how you expect. (I'm not sure what looks wrong -- sorry!) – Michael E2 May 18 '13 at 20:58
You can get the area contained by the contours with NIntegrate[Boole[PDF[dist, {x, y}] > 2.6*10^-6], {x, 0, 1024}, {y, 0, 1024}]. You might need to play with the options to get all the areas. – Michael E2 May 18 '13 at 21:02
Thanks I tried it out and this was the error message I got: NIntegrate::inumr: "The integrand Boole[PDF[dist,{x,y}]>1.152*10^-6] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1024},{0,1024}}" – Doron May 19 '13 at 8:50
Also, changing PlotRange parameter to "All" still didnt give the desired solution. (I wonder how I can put code lines and images in the comment area...) – Doron May 19 '13 at 9:01
Sorry, I realized used the variable name dist instead of \[ScriptCapitalD] for the distribution -- did you make the change back to script D? (You can edit the question to add more code/images that help explain the problem/question. That's probably better that putting them in a comment.) – Michael E2 May 19 '13 at 12:18

3 Answers

You can play with the method of integration. The default, "GlobalAdaptive", is not very good in this case. You should also consider how much precision is desired. For instance, the example distribution is based on a distribution of 3000 points and the resulting PDF is interpolated on 64 by 64 grid. I'll seek a precision of three digits.

"MultiPeriodic" works without any (helpful) warnings, but perhaps is not as accurate as requested (see further below):

NIntegrate[
Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024},
PrecisionGoal -> 3,
Method -> {"MultiPeriodic", "SymbolicProcessing" -> 0}] // Timing

{0.497470, 291727.}


Increase the requested precision to 4 and we get a warning:

NIntegrate[
Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024},
PrecisionGoal -> 4,
Method -> {"MultiPeriodic", "SymbolicProcessing" -> 0}] // Timing

NIntegrate::maxp: The integral failed to converge after 155092 integrand evaluations.
NIntegrate obtained 292514.50755723054 and 748.1766581075457 for the integral and
error estimates. >>

{2.205138, 292515.}


Note that the first answer is not within interval indicated by the second answer and error estimate.

With "LocalAdaptive", one can specify a partitioning of the domain. (No warning messages.)

NIntegrate[
Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024},
PrecisionGoal -> 3,
Method -> {"LocalAdaptive", "Partitioning" -> 128, "SymbolicProcessing" -> 0}] // Timing

{4.386049, 292851.}


Monte Carlo methods can work well in such cases, if high precision is not requried.

NIntegrate[
Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024},
Method -> {"QuasiMonteCarlo", "Partitioning" -> 64,
"MaxPoints" -> 10, "SymbolicProcessing" -> 0}] // Timing

{6.179461, 292700.}


Compare with SEngstrom's answer, which is basically a lower-endpoint Riemann sum:

(data2 = Table[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1023}, {y, 0, 1023}];
Select[Flatten[data2], # > 11.52*10^-7 &] // Length) // Timing

{8.264056, 292877}


All these give roughly the same answer (to three digits), the first one seeming a sort of outlier.

Another approach to the contour plot, building on Simon Wood's answer, is to use the interpolation grid itself for the plot-points grid. Unfortunately, I have to mask the plot range by hand, for the frame to look nice. :/

With[{ifunc =
First @ Cases[PDF[\[ScriptCapitalD], {x, y}], _InterpolatingFunction, -1, Heads -> True]},
Show[
ContourPlot[ifunc[x, y],
Evaluate[Sequence @@ MapThread[Prepend, {ifunc["Domain"], {x, y}}]],
PlotRange -> {{0, 1024}, {0, 1024}},
PlotRangePadding -> Scaled[0.02],
PlotPoints -> Most@Dimensions@ifunc["Grid"], MaxRecursion -> 3,
Contours -> {2.52*10^-6}],
Graphics[{White, Rectangle[{1024, -21}, {1045, 1045}],
Rectangle[{-21, 1024}, {1045, 1045}],
Rectangle[{-21, -21}, {0, 1045}],
Rectangle[{-21, -21}, {1045, 0}]}]
]
]


-

To force a numerical integration over the discrete Boole[] seems unsuitable to this problem. Why not sample the distribution back to a non-continuous representation and approximate the integral there? You can get fancier with the integration approximation if you need to.

data2 = Table[
Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 1024}];
Plus @@ Select[Flatten[data2], # > 11.52*10^-7 &]

(* 0.441366 *)


As pointed out by Michael E2 the requested property, the area, has a zeroth order approximation of

Select[Flatten[data2], # > 11.52*10^-7 &] // Length

-
This approximates the volume under the PDF (or equivalently, the probability), but the question asks for the area, which can be approximated by Length@Select[Flatten[data2], # > 11.52*10^-7 &] (== 292877) – Michael E2 May 20 '13 at 3:50
@MichaelE2 thanks for the correction. – SEngstrom May 20 '13 at 13:54

Here's one way to get a contour plot showing all the regions without going to very large numbers of plot points.

The PDF of the SmoothKernelDistribution is just an InterpolatingFunction, so we can examine the underlying data and extract the grid positions at which the function value exceeds the contour value. Those grid positions can then be fed into the contour plot as additional sample points.﻿

pts = With[{ifunc = First @
Cases[PDF[\[ScriptCapitalD], {x, y}], _InterpolatingFunction, -1, Heads -> True]},
Extract[ifunc["Grid"],
Position[ifunc["ValuesOnGrid"], z_ /; z >= 2.52*10^-6]]];

ContourPlot[
Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 1024},
PlotRange -> All, PlotPoints -> {50, pts}, Contours -> {2.52*10^-6}]


-
(+1) You beat me to it! I'd forgotten how to get the "ValuesOnGrid". – Michael E2 May 20 '13 at 18:03
No doubt about it you guys are just WIZARDS!!! Thanks a lot. Doron – Doron May 20 '13 at 19:17