Skip to main content

plotting - How: 2D scatterplots with quantitative density-dependent coloring


Consider the following scatterplot of a 50K-point dataset:


ListPlot[data,
AspectRatio -> Automatic, PlotRange -> {{x0, x1}, {y0, y1}},
ImageSize -> Small,
Frame -> True,
FrameTicksStyle -> Directive[FontOpacity -> 0, FontSize -> 0]]

Mathematica graphics


The color quickly saturates as one moves from the edge of the distribution to its center. As a result, most of the density information is lost.



One can remedy this slightly by assigning an opacity below 1 to the points:


ListPlot[data,
PlotStyle -> Opacity[0.05],
AspectRatio -> Automatic,
PlotRange -> {{x0, x1}, {y0, y1}}, ImageSize -> Small,
Frame -> True,
FrameTicksStyle -> Directive[FontOpacity -> 0, FontSize -> 0]]

Mathematica graphics


But this solution still has a couple of shortcomings:




  1. its dynamic range is still fairly narrow (even though it's wider than it was before); thus, most of the data cloud is still shows as saturated color;

  2. there's no explicit quantitative scale (e.g. a colorbar) tying colors (or in this case, shades) to densities;


The dynamic range problem could be solved by using more hues. This is what's routinely done when plotting flow cytometry data. For example:


enter image description here


(IMO, the plots in the last set would be improved if they included a color key, showing the correspondence between colors and densities.)


My question is how can I provide such quantitative density information in these scatterplots using Mathematica?



Answer



I think SmoothDensityHistogram (docs here) is what you are looking for:



data1 = RandomVariate[BinormalDistribution[{0, 0}, {2, 3}, 0.5], 100000];
data2 = RandomVariate[BinormalDistribution[{3, 4}, {2, 2}, .1], 100000];
data = data1~Join~data2;

This is just some random sample data. If you plot it using ListPlot, you obtain the "blob" you mentioned:


ListPlot[data, AspectRatio -> 1]

Mathematica graphics


Here is the same data presented with a smoothed 2D-histogram instead:


SmoothDensityHistogram[data, ColorFunction -> "TemperatureMap"]


Mathematica graphics




Data comparison: Jim Baldwin brought up a good point in comments regarding the need to compare multiple datasets, both visually and numerically. In that case, DensityHistogram may be the best bet. This function essentially is the discrete version of SmoothedDensityHistogram; the advantage in this context is the fact that it also has built-in tooltips whose value can be configured to report on distribution properties such as the total counts in each bin, probability, the value of the probability density function calculated from the data distributions, etc. In particular, this function may be most interesting because it can automatically generate legends for its data as shown below. Here is the documentation for DensityHistogram.


For instance, using the data above:


DensityHistogram[data, "Wand", "Count", 
ColorFunction -> "TemperatureMap",
ChartLegends -> Automatic
]


density histogram with tooltip and legend


Instead of "Count", one could also request the bin height to represent the PDF, CDF, etc. In this case I chose Wand binning among the built-in options because to me it seemed to offer the best compromise between fine-grained binning that reproduced the overall "shape" of the data, and execution time (ca. 7s on my machine). Knuth binning looked even better, but it took almost one minute to calculate on the same dataset!


In passing, I'd also like to mention that these *DensityHistogram functions seem to work very similarly under the hood, differing mostly in the way they present the data. In particular, my understanding is that both start by recovering a smooth kernel distribution from the existing data, using a Gaussian kernel by default.


Alternatively, other approaches focused on layering contour lines on top of a smooth density histogram have also been discussed in this question (Contour lines over SmoothDensityHistogram) to which Jim and others have contributed interesting answers.


Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...