# How to make a marginal distribution plot using DensityPlot?

How would you make a marginal distribution plot using DensityPlot? Adapting Sjoerd's answer from this question using Epilog, I can plot a slice at a given coordinate (for example, x=0 or y=0) as below:

f[x_, y_] := Exp[-2 (x^2 + y^2)] HermiteH[2, Sqrt x]^2;
DensityPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> All, Frame -> False,
Epilog -> {Line[Table[{x1, 2.05 + 0.2 f[x1, 0]}, {x1, -2, 2, 0.01}]],
Line[Table[{ 2.05 + 0.2 f[0, y1], y1}, {y1, -2, 2, 0.01}]] },
PlotRangePadding -> 0, PlotRangeClipping -> False, ImagePadding -> {{0, 100}, {0, 100}}]


This gives the following result: But what I'd really like is to plot the column- (row-) integrated values of DensityPlot along the x- (y-) axis margins.

The real function of interest for this calculation is expensive, so evaluating only once would be best (i.e., can we use the values of the DensityPlot?). Additionally, for this reason, DensityPlot is preferred over ListDensityPlot for its automatic mesh sampling, since the functions of interest tend to be localized, so a uniform mesh would be wasteful.

Below is an example where plotting along a given slice isn't representative of the marginal distribution: and the marginal plots (in red and blue, taken along the light dashed lines) just evaluate to zero.

• I think this is a duplicate of this question. – Jason B. Sep 10 '20 at 16:21
• As to your main problem, that the function f[x,y] is expensive to compute, and you don't want to have to keep recalling it for the density plot and the 1D plots, you have two choices. First, you could use memoization via f[x_, y_] := f[x,y] = <expensive computation>. Alternatively you could precompute the density as a matrix of values and call ListDensityPlot on it. – Jason B. Sep 10 '20 at 16:24
• @JasonB., great suggestion. ListDensityPlot is ~3x slower for similar resolution, which I chalked up to the adaptive mesh of DensityPlot, which I'd like to keep. Do you know how to extract the data (i.e., the (x,y,z) coordinates) underlying DensityPlot? If so, memoization would work really well in combination with @kglr 's answer below. – Jaffe42 Sep 10 '20 at 17:41
• do you have access to the chat room? – Jason B. Sep 10 '20 at 19:09

ClearAll[f, xMargin, yMargin, ppX, ppY]

f[x_, y_] := Exp[-2 (x^2 + y^2)] HermiteH[2, Sqrt x]^2

xMargin[x_] = Integrate[f[x, y], {y, -Infinity, Infinity}];

yMargin[y_] = Integrate[f[x, y], {x, -Infinity, Infinity}];

xrange = {-3, 3};

yrange = {-2, 2};

scale = 1/4/Pi;

gap = 0.05;

dp = DensityPlot[f[x, y], {x, xrange[], xrange[]}, {y, yrange[], yrange[]},
PlotRange -> All]


We can construct appropriately translated margins using ParametricPlot:

ppY = ParametricPlot[{xrange[] - gap - scale v yMargin[y], y},
{y, yrange[], yrange[]}, {v, 0, 1},
PlotStyle -> Red, PlotPoints -> 50, Axes -> False];

ppX = ParametricPlot[{x, yrange[] - gap - scale v xMargin[x] },
{x, xrange[], xrange[]}, {v, 0, 1},
PlotStyle -> Blue, PlotPoints -> 50, Axes -> False];


and combine them with dp using Show:

Show[ppY, ppX, dp, PlotRange -> All, Frame -> True] To show the marginal plots on top and right frames:

ppY2 = ParametricPlot[{xrange[] + gap + scale v yMargin[y], y},
{y, yrange[], yrange[]}, {v, 0, 1},
PlotStyle -> Red, PlotPoints -> 50, Axes -> False];

ppX2 = ParametricPlot[{x, yrange[] + gap + scale v xMargin[x]},
{x, xrange[], xrange[]}, {v, 0, 1},
PlotStyle -> Blue, PlotPoints -> 50, Axes -> False];

Show[ppY2, ppX2, dp, PlotRange -> All, Frame -> True] To put the marginal plots outside the frame, we can use Inset + Epilog:

insetY = Inset[#, {xrange[] (1 + gap), yrange[]},
{Left, Top}, Scaled] & @ ppY2;

insetX = Inset[#, {xrange[], yrange[] (1 + gap)},
{Right, Bottom}, Scaled] & @ ppX2;

Show[dp, Epilog -> {insetX, insetY},
ImagePadding -> {{Scaled[.02], Scaled[.1]}, {Scaled[.02], Scaled[.1]}},
ImageSize -> Large, PlotRangeClipping -> False, ] Alternatively, we can Plot the functions xMargin and yMargin and use GeometricTransformation with appropriate transformation functions position them and Show the transformed graphics objects with dp:

ClearAll[transform, tFX, tFY]
transform[tf_] := Graphics[#[] /.
ll : (_Line | _Polygon) :> GeometricTransformation[ll, tf]] &;

tFY = TranslationTransform[{-gap, xrange[]}]@*
RotationTransform[Pi/2, {xrange[], 0}];

tFX = TranslationTransform[{0, yrange[] - gap}]@*
ScalingTransform[{1, -1}];

pltY = Plot[scale yMargin[y], {y, yrange[], yrange[]},
Filling -> Axis, PlotStyle -> Red, Axes -> False];
pltX = Plot[scale  xMargin[x], {x, xrange[], xrange[]},
Filling -> Axis, PlotStyle -> Blue, Axes -> False];

Show[transform[tFY]@pltY, transform[tFX]@pltX, dp, PlotRange -> All,
Frame -> True] To show the marginal plots on top and right frames use the transformations tFX2 and tFY2:

tFY2 = TranslationTransform[{gap, xrange[]}]@*
RotationTransform[-Pi/2, {xrange[], 0}];

tFX2 = TranslationTransform[{0, yrange[] + gap}];

Show[transform[tFY2] @ pltY, transform[tFX2] @ pltX, dp, PlotRange -> All,
Frame -> True] Update: An alternative approach to get the marginal plots: Use Plot3D to plot f with equally spaced mesh lines in x and y directions and extract the coordinates of mesh lines.

ndivs = 50;
{meshx, meshy} = Subdivide[#[], #[], ndivs] & /@ {xrange, yrange};

coords = Plot3D[f[x, y],
{x, xrange[], xrange[]}, {y, yrange[],  yrange[]},
PlotRange -> All, Mesh -> {meshx, meshy}, PlotStyle -> None][[1, 1]];


Group coords by the first and second coordinates and construct two WeightedData objects and plot them using SnoothHistogram:

bw = .01;
{wDx, wDy} = Table[Apply[WeightedData] @ Transpose @ KeyValueMap[List] @
KeySort @ GroupBy[coords, Round[#[[i]], bw] & -> Last, Mean], {i, 2}];

{sHx, sHy} = {SmoothHistogram[wDx, PlotStyle -> Blue,
Filling -> Axis, ImageSize -> 300],
SmoothHistogram[wDy, PlotStyle -> Red, Filling -> Axis, ImageSize -> 300]};

Row[{sHx, sHy}, Spacer] Alternatively, Plot the PDF of SmoothKernelDistribution of wDx and wDy:

{sKDx, sKDy} = SmoothKernelDistribution /@ {wDx, wDy};

{sHx2, sHy2} = {Plot[PDF[sKDx]@x, {x, xrange[], xrange[]},
PlotStyle -> Blue, Filling -> Axis, ImageSize -> 300],
Plot[PDF[sKDy]@y, {y, xrange[], yrange[]}, PlotStyle -> Red,
Filling -> Axis, ImageSize -> 300]};

Row[{sHx2, sHy2}, Spacer] Update 2: Processing DensityPlot output to get {x,y,z} coordinates (where z is scaled to the unit interval:

dp = DensityPlot[f[x, y], {x, -3, 3}, {y, -2, 2},
ColorFunction -> Hue, PlotRange -> All, PlotPoints -> 50]
coordsFromDP = Join[dp[[1, 1]], List /@ dp[[1, 3, 2, All, 1]], 2];


Except for the scale of the z coordinate ListPlot3D of coordsFromDP is "close" to the Plot3D output:

Row @ {Plot3D[f[x, y], {x, -3, 3}, {y, -2, 2}, ImageSize -> 300,
PlotRange -> All], ListPlot3D[coordsFromDP, ImageSize -> 300]} We process coordsFromDP the same way we did for coords above (except for a larger bin width):

bw = .02;
{wDx2, wDy2} = Table[Apply[WeightedData] @ Transpose @ KeyValueMap[List] @
KeySort@GroupBy[coordsFromDP, Round[#[[i]], bw] & -> Last,  Mean], {i, 2}];

{sHx2, sHy2} = {SmoothHistogram[wDx2, PlotStyle -> Blue,
Filling -> Axis, ImageSize -> 300],
SmoothHistogram[wDy2, PlotStyle -> Red, Filling -> Axis, ImageSize -> 300]};

Row[{sHx2, sHy2}, Spacer] • Wow, this is really beautiful! Unfortunately, the real function I'm interested takes quite a long time to evaluate. Do you know a way of using the data within dp to generate pltX and pltY? If so, that would plug right into your GeometricTransformation positioning/formatting approach and fully solve my problem!! – Jaffe42 Sep 10 '20 at 14:32
• @Jaffe42, I updated with an alternative approach using Plot3D output to generate the marginal plots. – kglr Sep 10 '20 at 16:09
• Wow, this is amazingly comprehensive. Completely solved the problem from all angles. Thank you so much! – Jaffe42 Sep 10 '20 at 20:31