Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 5:32
1  
Please include an example in your question. Do not leave us guessing. –  Mr.Wizard Mar 14 '12 at 5:35
    
Also, welcome to Mathematica.SE. –  Mr.Wizard Mar 14 '12 at 5:36
    
Thanks Mr. Wizard. –  Henry B Mar 14 '12 at 6:14
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 7:04
show 4 more comments

2 Answers

up vote 8 down vote accepted

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

share|improve this answer
    
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
show 3 more comments

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)

share|improve this answer
    
+1 The last plot is actually a reflection, not a rotation. –  whuber Mar 16 '12 at 14:01
    
@whuber True, a reflection in the x==y axis. –  Sjoerd C. de Vries Mar 16 '12 at 14:57
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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