A function to plot the confidence ellipse of the covariance of a 2D dataset. Uses matplotlib.

plot_confidence_ellipse.py

import numpy as np
import matplotlib.pyplot as plt
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`
    
    See how and why this works: https://carstenschelp.github.io/2018/09/14/Plot_Confidence_Ellipse_001.html
    
    This function has made it into the matplotlib examples collection:
    https://matplotlib.org/devdocs/gallery/statistics/confidence_ellipse.html#sphx-glr-gallery-statistics-confidence-ellipse-py
    
    Or, once matplotlib 3.1 has been released:
    https://matplotlib.org/gallery/index.html#statistics
    
    I update this gist according to the version there, because thanks to the matplotlib community
    the code has improved quite a bit.

    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.

    Returns
    -------
    matplotlib.patches.Ellipse

    Other parameters
    ----------------
    kwargs : `~matplotlib.patches.Patch` properties
    """
    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-dimensionl 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 stdandard 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 stdandard 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)
    # render plot with "plt.show()".

[...] Thank you again! This was exactly what I needed!

Hi @ktevans, This is great. You turned it into your own thing.
You didn't even need the arctan(), after all. The slope was all you needed.
With the sign fixed, then.
I am pleased to see that this was useful.
Kind regards, Carsten

Hello Carsten,
Firstly, thank you for sharing this script and for all the information concerning it. I was wondering if it is possible to extract the angle of the ellipse (either in radians or degrees). I have already the radii and the width (minor and major). It would be awesome if I was able to get the "tilt" / the angle of the ellipse!
Thanks a million.
Best regards, Antoine Dupuy.

This is also exactly what I am trying to get! Still working on it, but if anyone has any insight, I'd appreciate it!

Hi @ktevans and hi @Antoine-DupuyUL,

(And sorry Antoine, your question from months ago must have gotten lost in the shuffle). I am trying something from top of my head, here. I am not at home with not enough time.

In a normalized case, that is, when the standard deviations of both variables are the same, the angle of the ellipse will be

arctan(std1/std2) = arctan(1) = pi/4 = 45degrees.

In the article I also mention that any deviation from 45degrees only comes forth from different scaling of the variables. So you might give it a shot to determine the angle by

arctan(std1/std2) = angle in rad

std1 and std2 would be the "original", non-normalized standard deviations. This will alsways be an "upward" angle. A seemingly "downward" angle that you see when the pearson coefficient is negative is only due to the "second" radius being longer than the "first" radius. Does this work for you?

Kind regards, Carsten

Hello Carsten,

Thanks a million for your help and all the explanation!
It works perfectly 🙏

Again, thank you for making this code available.

Best regards,
Antoine.