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()".

If you (and others) agree, I can make a PR over at Matplotlib. Cheers.

Hi @liamtoney
You are more than welcome to fix that on the Matplotlib page!
Actually, I made an effort to straighten out this flaw, end of last year.
But the review process turned out to be less supportive than I had experienced it earlier.
(Generally, the Matplotlib crew is very positive and enthusiast!)
The agreement that I remember was that the maintainers themselves would submit a minimal change that would correct the wrong statement.
But it seems like it didn't happen.
So if you independently make the same point and can point to an authoritative source, it should be no problem to submit a PR that passes smoothly.
Thanks a lot!
Regards, Carsten

Hi @adiringenbach,
Even this short answer took me a long time. Sorry about that.
I think that a lot of the confusion comes from the fact that I used the word "vertices" to notify only the extreme points of an ellipse.
Whereas in graphicsland "vertex" is just another word for point. More or less.
I have yet an alternative method for getting (the points of) a covariance ellipse.
The documentation is not as detailed as in this gist, though:

In this new gist I am sketching the principle.
I also explain how it can be combined with the code in this current gist and become a complete and usable function.
Hope this is still of any use to you or someone else.

Here's the merged PR with the correct confidence interval: matplotlib/matplotlib#21216

Hi Carsten,
thanks for implementing this! I am using your code in a script and I'd like to cite your work properly. AFAIK Github now offers a way to add a "cite" button to repos (?) as well as getting a DOI. Would you consider getting one set up for this code? Otherwise, I'd just cite the URL, but I always prefer somwthing official as it offers a more consistent solution...
Thanks,
Jens

Hi @jensweiser,
Thank you. DOI sounds interesting, anyway. I will see what I can do.

Hi @jensweiser,

That was quite easy.
Obviously, I cannot get a DOI for a gist, but the repo that feeds the github.io page has a release with a zenodo DOI, now.
The file you want to refer to would be _posts/2018-09-14-Plot_Confidence_Ellipse_001.md
Here's the data of the DOI:

Markdown:
[![DOI](https://zenodo.org/badge/138745361.svg)](https://zenodo.org/badge/latestdoi/138745361)

reStructedText:

.. image:: https://zenodo.org/badge/138745361.svg
   :target: https://zenodo.org/badge/latestdoi/138745361

or as HTML:
<a href="https://zenodo.org/badge/latestdoi/138745361"><img src="https://zenodo.org/badge/138745361.svg" alt="DOI"></a>

Hi @jensweiser,

That was quite easy. Obviously, I cannot get a DOI for a gist, but the repo that feeds the github.io page has a release with a zenodo DOI, now. The file you want to refer to would be _posts/2018-09-14-Plot_Confidence_Ellipse_001.md Here's the data of the DOI:

Markdown: [![DOI](https://zenodo.org/badge/138745361.svg)](https://zenodo.org/badge/latestdoi/138745361)

reStructedText:

.. image:: https://zenodo.org/badge/138745361.svg
   :target: https://zenodo.org/badge/latestdoi/138745361

or as HTML: <a href="https://zenodo.org/badge/latestdoi/138745361"><img src="https://zenodo.org/badge/138745361.svg" alt="DOI"></a>

....fantastic, thanks a lot!

Hi SpartyOn, Thank you for your question. To be honest I have no idea what exactly the returned object of the function contains. As I mentioned earlier in these comments, I was advised to return the ellipse object because it might be useful for some users. I have no idea if the height and width of it are of any use, here. To get the actual radii of the ellipse at scale you'd have to write some code that parallels the "transforms" part of the function. That means, emulates the rotating and scaling. Not elegant but it will give you that "a" and "b" of the ellipse. I have created another gist with a modified function that also returns the radii of the ellipse. Its this gist: https://gist.github.com/CarstenSchelp/8de3924c7008eb3f5894549f0f04cb96 As I said, not pretty, but might work. I just tested it "by sight" no guarantees given whatsoever! But you can see the approach to get the values and can do some debugging if necessary. Hope that this helps. Kind regards, Carsten

…

On 4/6/22 21:56, SpartyOn wrote: ***@***.**** commented on this gist. ------------------------------------------------------------------------ How would you go about calculating the area of the ellipse you make? I tried this, but it's giving me values that are way too small. |ellipse = confidence_ellipse(df['x'], df['y'], ax, alpha=0.5, facecolor='green', edgecolor='green', zorder=0) ax.scatter(df['x'], df_test['y'], s=0.5,c='brown') ell_height_radius = (ellipse.get_height())/2 ell_width_radius = (ellipse.get_width())/2 ellipseArea = np.pi*ell_height_radius*ell_width_radius plt.show()| — Reply to this email directly, view it on GitHub <https://gist.github.com/b992645537660bda692f218b562d0712#gistcomment-4123919>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD67ZVCHNFBRPTYANFDRJJTVDXT5TANCNFSM4ITJFVGA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

Do you think is it possible to apply a variant of this to 3D coordinates (xyz)?

Hi!
I tested the gist to get the values of the radii and it doesn't seem to give the good values. I spent a bit of time trying to do some debugging but I can't figure out the problem! Did anyone have the same problem?

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.

Hi! I tested the gist to get the values of the radii and it doesn't seem to give the good values. I spent a bit of time trying to do some debugging but I can't figure out the problem! Did anyone have the same problem?

If you are taking the "ell_radius_x" & "ell_radius_x" directly, you will find weird results. That is because they are not scaled to your data yet.
Try with this:
scaled_ell_x_radius = np.sqrt(cov[0, 0]) * n_std
scaled_ell_y_radius = np.sqrt(cov[1, 1]) * n_std

There should be more logical results!

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!

Do you think is it possible to apply a variant of this to 3D coordinates (xyz)?

Guess what I have been banging my head against, the last six years :D
It's tricky, to say the least. 2D is just such a nice special case.
Don't hold your breath, but stubborn as I am, I keep trying.
Which is bold, because I think that people with profoundly deeper math-chops than I have tried this, too.
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,
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

Hi Carsten,

Thank you for this code. It has been of great help to me! I was able to get the eccentricity of the ellipse using this method!

I defined a simple correlation variable via the following:

corr = 0
if pearson > 0:
  corr = 1
if pearson < 0:
  corr = -1

Then, the eccentricity is just:
ecc = corr * (np.std(y) / np.std(x))

I got the following results where the slope in the legend is the calculated eccentricity:

Thank you again! This was exactly what I needed!

[...] 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.