# How to plot a multivariate normal distribution?

**URL:** <https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486>\
**Category:** Visualization\
**Created:** [March 3, 2020, 10:04pm UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486 "2020-03-03T22:04:03Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [March 3, 2020, 10:04pm UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/1 "2020-03-03T22:04:03Z")

</div>

Plotting univariate distributions is really simple:

```julia
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:

```julia
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:

 ![](https://global.discourse-cdn.com/julialang/original/3X/2/8/281c590d6bcf211857306a0fbd44634f9f38c7bc.png)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [March 3, 2020, 10:17pm UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/2 "2020-03-03T22:17:16Z")

</div>

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.

```julia
Z = [pdf(mvnorm,[i,j]) for i in 0:100, j in 0:100]
plot(0:100,0:100,Z,st=:surface)

```

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [March 3, 2020, 10:23pm UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/3 "2020-03-03T22:23:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [March 4, 2020, 1:18am UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/4 "2020-03-04T01:18:53Z")

</div>

> [@mthelm85](#):
>
> Z = [pdf(mvnorm,[i,j]) for i in 0:100, j in 0:100]  
> plot(0:100,0:100,Z,st=:surface)

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! 😆 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.

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [March 4, 2020, 1:45am UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/5 "2020-03-04T01:45:10Z")

</div>

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/1/b/1b25596aa11daeb55dc71b8ef2a548a2fcec0da3.png)

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [March 4, 2020, 6:50am UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/6 "2020-03-04T06:50:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Robert\_Moss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robert_moss/32/48004_2.png) [@Robert\_Moss](https://discourse.julialang.org/u/Robert_Moss)\
**Post date:** [February 23, 2021, 11:37pm UTC](https://discourse.julialang.org/t/how-to-plot-a-multivariate-normal-distribution/35486/7 "2021-02-23T23:37:34Z")

</div>

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:

```julia
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:

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

```

![image](https://global.discourse-cdn.com/julialang/original/3X/b/0/b067d60d9b8474c9a1df3e28aa6b9c4a8282e48a.png)