# How to make a discrete density plot for a two-variable PDF?

I have the following function of two variables

$$(1)\qquad f(x,y)=\exp[-x^{2} - y^{2} + \frac{1}{2}xy]$$

which is interpreted as a two-variable Gaussian probability distribution. Concretely I want to make a discrete density plot in the $$xy$$-plane for the probability function (1), instead of the continuous ones usually constructed by the DensityPlot.

Is there any way to make a discrete density plot for the probability distribution (1)? I mean a set of points that are distributed according to the probability distribution (1).

• I don't understand what you mean by a "discrete" density plot. – mikado Dec 15 '19 at 21:20
• @mikado By a "discrete" density plot, I mean a set of points that are distributed according to the probability distribution (1). – Julio Abraham Mendoza Fierro Dec 15 '19 at 21:32
• Would you explain why you'd want a discrete representation of a smooth surface? Also, an answer to a question posted by others: Do you have a set of samples or do you have a function describing a surface? Both answers would be helpful. – JimB Dec 16 '19 at 16:28

Is this what you mean? When you "a set of points distributed according to..." it's not clear if you want them randomly sampled from that or what.

dist = MultinormalDistribution[{0, 0}, 1/3*{{2, 1}, {1, 2}}];

PDF[dist, {x, y}] // Simplify


$$\frac{\sqrt{3} e^{-x^2+x y-y^2}}{2 \pi }$$

Histogram3D[RandomVariate[dist, 10000]]


Clear["Global*"]

distbn = BinormalDistribution[{0, 0}, {s1, s2}, p];

PDF[distbn, {x, y}] // Expand

(* E^(-(x^2/(2 (1 - p^2) s1^2)) + (p x y)/((1 - p^2) s1 s2) - y^2/(
2 (1 - p^2) s2^2))/(2 Sqrt[1 - p^2] π s1 s2) *)


Solving for {s1, s2, p} to match the specified form

Solve[{2 (1 - p^2) s1^2 == 1, 2 (1 - p^2) s2^2 == 1, ((1 - p^2) s1 s2)/p == 2,
s1 > 0, s2 > 0}, {s1, s2, p}]

(* {{s1 -> 2 Sqrt[2/15], s2 -> 2 Sqrt[2/15], p -> 1/4}} *)


The distribution is

dist = BinormalDistribution[{0, 0}, {2 Sqrt[2/15], 2 Sqrt[2/15]}, 1/4];

PDF[dist, {x, y}] // Simplify

(* (Sqrt[15] E^(-x^2 + (x y)/2 - y^2))/(4 π) *)

SeedRandom[1234]

data = RandomVariate[dist, 2000];

Needs["MultivariateStatistics"]

With[{q = {0.5, 0.9, 0.95}},
Show[
ListPlot[data, PlotStyle -> Red],
Graphics /@
EllipsoidQuantile[data, q]]]