How to plot a multivariate normal distribution?

Domains Visualization
#1

Plotting univariate distributions is really simple:

using Distributions
using StatsPlots

norm = fit(Normal, rand(100))
plot(norm)

However, how would I plot a multivariate distribution? The above code doesn’t work for the multivariate case:

mvnorm = fit(MvNormal, [rand(0.0:100.0, 100) rand(0.0:100.0, 100)]')

julia> plot(mvnorm)
ERROR: MethodError: no method matching iterate(::MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}})
Closest candidates are:
  iterate(::Core.SimpleVector) at essentials.jl:600
  iterate(::Core.SimpleVector, ::Any) at essentials.jl:600
  iterate(::ExponentialBackOff) at error.jl:218
  ...

Here’s an example from Wikipedia of what an ideal output would look like:

1280×800 35.9 KB

2 Likes
#2

Doesn’t look like there are any shortcuts for it in StatsPlots.

You can just get the pdf values on a grid and do a surface plot though.

Z = [pdf(mvnorm,[i,j]) for i in 0:100, j in 0:100]
plot(0:100,0:100,Z,st=:surface)
7 Likes
#3

Yes. It’d be nice to add a recipe for that to StatsPlots

2 Likes
#4

This is great, thanks! Any idea how to change the color of the grid lines when doing st=:wireframe? I’ve tried every single color keyword argument/alias I can find in the docs and none of them work! :laughing: Or, when using st=:surface do you know how to change the alpha? Again, I’ve tried all the alpha-related kwargs and none seem to do the trick.

#5

I was able to hack together something that looks decent. Below is the code as well as the output for some data I’m working with:

x = 0:maximum(X)
y = 0:(maximum(Y)/length(x))+1:maximum(Y) # because x and y are of different lengths
z = [pdf(mvnorm, [i, j]) for i in x, j in y]

pyplot()

plot(x, y, z, linetype=:surface, legend=false, color=:blues)
plot!(x, y, z, linetype=:wireframe, legend=false, color=:black, width=0.1)

test

2 Likes
#6

Looks nice. wireframe color should follow the linecolor keyword. alpha for surface isn’t implemented :disappointed:
BTW no need to call collect all the time (and you can delete the first assignment to y).

1 Like
#7

Just a warning/note about x-y ordering when creating the Z grid: if you’re plotting as a contour/heatmap, then you’ll want to reverse the for x in X, y in Y ordering to be for y in Y, x in X so that the grid maps to the x-y axes correctly. This is subtle, as it will not always present itself as a problem depending on the parameters of the MvNormal.

Here’s an example plotting a bivariate normal distribution as filled contours:

using Distributions
using Plots

μ = [0, 0]
Σ = [1  0.9;
     0.9 10]
p = MvNormal(μ, Σ)

X = range(-8, 8, length=100)
Y = range(-8, 8, length=100)
Z = [pdf(p, [x,y]) for y in Y, x in X] # Note x-y "for" ordering
contourf(X, Y, Z, color=:viridis)

To avoid Z altogether, you can pass a function f to handle calculating the z-values:

f(x,y) = pdf(p, [x,y])
contourf(X, Y, f, color=:viridis)

image

2 Likes