Hi Matplotlib folks,
I’m trying to add contour lines to a roughly bivariate normal scatter plot at fixed standard deviations away. There’s a very nice page in the documentation on how to do this by finding the right translation and transformation and drawing an ellipse, but this does not allow for the same flexibility in labeling as ax.contour() and ax.clabel(). But ax.contour() also is an odd solution as well, since it’s interpolated from a grid, and the exact form of the ellipse is known. Before I embark on trying to figure out how to do this from scratch: what is easier? Getting ax.contour() to give me the standard error contours, and label with ax.clabel()? Or add labels to the ellipse somehow?
Below is some code showing the method from the documentation
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import multivariate_normal as mvn
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms
def confidence_ellipse(x, y, ax, n_std=3.0, facecolor='none', **kwargs):
"""
Create a plot of the covariance confidence ellipse of *x* and *y*.
Parameters
----------
x, y : array-like, shape (n, )
Input data.
ax : matplotlib.axes.Axes
The axes object to draw the ellipse into.
n_std : float
The number of standard deviations to determine the ellipse's radiuses.
**kwargs
Forwarded to `~matplotlib.patches.Ellipse`
Returns
-------
matplotlib.patches.Ellipse
"""
if x.size != y.size:
raise ValueError("x and y must be the same size")
cov = np.cov(x, y)
pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
# Using a special case to obtain the eigenvalues of this
# two-dimensional dataset.
ell_radius_x = np.sqrt(1 + pearson)
ell_radius_y = np.sqrt(1 - pearson)
ellipse = Ellipse((0, 0), width=ell_radius_x * 2, height=ell_radius_y * 2,
facecolor=facecolor, **kwargs)
# Calculating the standard deviation of x from
# the squareroot of the variance and multiplying
# with the given number of standard deviations.
scale_x = np.sqrt(cov[0, 0]) * n_std
mean_x = np.mean(x)
# calculating the standard deviation of y ...
scale_y = np.sqrt(cov[1, 1]) * n_std
mean_y = np.mean(y)
transf = transforms.Affine2D() \
.rotate_deg(45) \
.scale(scale_x, scale_y) \
.translate(mean_x, mean_y)
ellipse.set_transform(transf + ax.transData)
return ax.add_patch(ellipse)
mean, cov = (2, 1), ((1, 0.5), (0.5, 1))
dat = mvn(mean, cov).rvs(1000)
fig, ax = plt.subplots()
ax.scatter(*dat.T)
for i in range(1, 3):
confidence_ellipse(*dat.T, ax, edgecolor='r', n_std=i)
thank you,
Vince