10

How do I acheive something like InterpolationOrder in ListContourPlot for ContourPlot? My main problem is I want to smooth the contours in ContourPlot.

I have some points which I'm going to call data, which is 10,000 points and looks like this:

{{24.3,40.0},{29.2,56.0} ... }

I create this distance function:

D[x_,y_]=(#[[1]]-x)^2+(#[[2]]-y)^2&/@data

Now, my function is supposed to count the number of points within a certain distance of any x,y position and divide by the total number of points:

F[x_,y_,r_]:=Length[Select[D[x,y],(#<r^2)&]]/Length[data]

Now I can plot my function,

ContourPlot[F[x,y,5],{x,0,100},{y,0,100}]

The problem is my function is not naturally smooth.

I know about SmoothDensityHistogram and DensityHistogram, but I like the more simple look of smoothed contours. Can anyone help with this?

Sjoerd C. de Vries
  • 65,815
  • 14
  • 188
  • 323
Henry B
  • 447
  • 1
  • 3
  • 12
  • The problem is ContourPlot takes a function. I would use ListContourPlot, but I want to superimpose another plot on top of the contour and I lose the (x,y) information with ListContourPlot. – Henry B Mar 14 '12 at 05:32
  • 1
    Please include an example in your question. Do not leave us guessing. – Mr.Wizard Mar 14 '12 at 05:35
  • Also, welcome to Mathematica.SE. – Mr.Wizard Mar 14 '12 at 05:36
  • Is there a reason you aren't using EuclideanDistance[#,{x,y}]&/@ data and are instead rolling your own distance function? – Verbeia Mar 14 '12 at 06:35
  • @Verbela: Yes. The reason is I'm not very good with Mathematica. :) Thanks for the suggestion. I will adjust my code. – Henry B Mar 14 '12 at 06:43
  • 4
    One general suggestion: Don't define functions/symbols that start with a capital letter as they might shadow or clash with MMa defined functions and symbols. An example is your use of D which MMa is using for the derivative. – Matariki Mar 14 '12 at 07:04
  • @szabolcs the function Henry wants to plot is a discrete one, but he wants the steps to be smoothed. Hence, the need for some kind of interpolation. – Sjoerd C. de Vries Mar 14 '12 at 07:30
  • @Szabolcs Ah, I see. I didn't go that far back in the edits. – Sjoerd C. de Vries Mar 14 '12 at 08:01

2 Answers2

8

The most efficient way to carry out your task, which is to plot a contour map of a kernel density of your points, is by converting the points to raster format and using a Fast Fourier Transform to convolve them with a density kernel. But that takes some work. If you're willing to wait a few seconds, the whole procedure is (less efficiently) built into Mathematica's SmoothKernelDistribution function.

Here is an example taken, with minor changes (to make it more interesting), directly from the help page:

(* Create some data--around 10,000 points--for the illustration *)
data = Join @@ Table[RandomVariate[BinormalDistribution[m, {1/2, 1/2}, 0], 1500], 
        {m, RandomReal[{1, 9}, {7, 2}]}];

(* Create a rough (D1) and smooth (D2) density for contouring *)
D1 = SmoothKernelDistribution[data, 0.02]; (* Takes a few seconds *)
D2 = SmoothKernelDistribution[data, 0.5];  (* Takes a few more seconds *)

(* Plot the points and their densities *)
points = ListPlot[data, PlotRange -> {{0, 10}, {0, 10}}, AspectRatio -> 1];
TableForm[{
  Prepend[
   Table[
     ContourPlot[
       Evaluate@PDF[D, {x, y}], {x, 0, 10}, {y, 0, 10}, 
       PlotRange -> All, 
       ColorFunction -> "TemperatureMap"], 
     {D, {D1, D2}}
 ], points]
}]

Figure

whuber
  • 20,544
  • 2
  • 59
  • 111
  • Actually your code runs near instantaneously here. Did you use a lot more points when you found that it takes a few seconds? If you don't want to use SmoothKernelDistribution because of performance concerns (actually I don't know how it's implemented...), and you'd rather use something FFT based, that's easy to do as well: ListConvolve[GaussianMatrix[5], BinCounts[data, .1, .1]] // ListContourPlot – Szabolcs Mar 14 '12 at 16:58
  • Thanks. I don't think ListConvolve, as directly applied to the data, does what the OP intends. You first need to create a 2D array, somewhat as in your reply here, but you don't evaluate the density F on a regular grid: you grid the points (as sums of indicator functions) and then you could use ListConvolve et al. with a suitably chosen kernel. SmoothKernelDistribution really bogged down for me when using custom kernels; at one point it ate all my RAM and forced a reboot :-(. – whuber Mar 14 '12 at 17:02
  • I see you're using BinCounts to rasterize the points, so that's ok: nice, simple solution. In practice, far more bins would likely be desirable. – whuber Mar 14 '12 at 17:04
  • "The most efficient way to carry out your task, which is to plot a contour map of a kernel density of your points, is by converting the points to raster format and using a Fast Fourier Transform to convolve them with a density kernel." <-- Doesn't this require a regular binning? Is there a way to use the FFT to speed up the convolution when we don't have a function sampled on a regular grid? Yes, you are right that the BinCounts + ListConvolve method will get slow quickly as the bin size is decreased. – Szabolcs Mar 14 '12 at 17:24
  • It might be worth some testing and benchmarking. In principle, using sparse matrices, BinCounts is just a $O(n)$ operation (time and storage) for $n$ input points. ListConvolve--or at least the FFT equivalent--will need $O(1/h^2)$ storage and $O(-\log(h)/h)$ time for a bin width of $h$. – whuber Mar 14 '12 at 17:46
  • I assume smoothing influences the height levels of the plot, for instance, its extreme values, right? So, wouldn't that be an unwanted side effect of creating smoothed contours? – Sjoerd C. de Vries Mar 16 '12 at 14:56
  • Yes, it is worth emphasizing that smoothing affects height--otherwise it couldn't possibly change the contours! Whether this unavoidable effect is "unwanted" depends on whether one is aware of what smoothing does. – whuber Mar 16 '12 at 16:09
  • @whuber out of curiosity do you recall what custom kernels caused the system to "bog down"? – Andy Ross Apr 14 '12 at 05:41
6

I'd say you could solve this problem by evaluating your function F on a regular grid and feeding this data to ListContourPlot with the InterpolationOrder set at an appropriate level.

data = RandomVariate[
          MultinormalDistribution[{50, 50}, 150 {{1, 0.8}, {0.8, 4}}], 10000];
dist[x_, y_] := EuclideanDistance[#, {x, y}] & /@ data;
f[x_, y_, r_] := Length[Select[dist[x, y], (# < r^2) &]]/Length[data]
pts = Table[f[x, y, 5], {x, 0, 100}, {y, 0, 100}]; 

ListContourPlot[pts, InterpolationOrder -> 0]

Mathematica graphics

ListContourPlot[pts, InterpolationOrder -> 1]

Mathematica graphics

The job ContourPlot does isn't that bad, by the way:

ContourPlot[f[x, y, 5], {x, 0, 100}, {y, 0, 100}]

Mathematica graphics

(Yes, I know the plot is rotated with respect to the above ones. I forgot to order x and y in the correct order in the Table function)

Sjoerd C. de Vries
  • 65,815
  • 14
  • 188
  • 323