Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated scaling factors for proper geometric representation of uncertainty ellipses #1067

Merged
merged 4 commits into from
Apr 12, 2022
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 27 additions & 15 deletions python/gtsam/utils/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,9 @@ def plot_covariance_ellipse_3d(axes,
Plots a Gaussian as an uncertainty ellipse

Based on Maybeck Vol 1, page 366
k=2.296 corresponds to 1 std, 68.26% of all probability
k=11.82 corresponds to 3 std, 99.74% of all probability
For the 3D case:
k = 3.527 corresponds to 1 std, 68.26% of all probability
k = 14.157 corresponds to 3 std, 99.74% of all probability
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have an online reference?

Copy link
Contributor

@senselessDev senselessDev Jan 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With code from other comment:

pct_to_sigma(sigma_to_pct(1, 1), 3) ** 2 --> 3.5267403802617303
pct_to_sigma(sigma_to_pct(3, 1), 3) ** 2 --> 14.156413609126677

So seems about right if we want this behavior.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these (3.5 and 14.1) should still be changed to standard deviation


Args:
axes (matplotlib.axes.Axes): Matplotlib axes.
Expand All @@ -87,7 +88,7 @@ def plot_covariance_ellipse_3d(axes,
n: Defines the granularity of the ellipse. Higher values indicate finer ellipses.
alpha: Transparency value for the plotted surface in the range [0, 1].
"""
k = 11.82
k = np.sqrt(14.157)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

above = 14.157 and here it is sqrt(14.157)

Did you take a look at the PR #1063 as I suggested?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I took a look at PR #1063 but am not sure how it is relevant. Perhaps I am missing something? I am new to contributing to open source projects so perhaps there was a file or something I needed to look at. As far as I can tell, the discussion there is regarding removing a negative sign in equations 5.3 and 5.6 in the math.pdf document.

As for the k values I chose, they come from the Chi Squared distribution's p-table for different dimensions. E.g., for a single dimension, k = 9 corresponds to 3 standard deviations, for 2D we have k = 11.820, for 3D we have k = 14.157, and for 4D k = 16.251 (all values taken from the table on page 366 of the Maybeck Vol 1 textbook). Note that I am choosing k to be the variance rather than the standard deviation, which is why you see the square root. I chose to write the comments as I did simply because I wanted to keep them as similar as possible to the previous comments pointing out the k values.

The most useful online reference I could find for choosing these k values can be found here and the Maybeck textbook can be found here.

Please let me know if I can answer any other questions.

Thanks!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Take a look at my note at the top of plot.py. I added some python code for 2D.

  1. In 2D I chose to leave the k=5, because it is consistent with what you said in the comments, and it is referenced in the documentation.
  2. I don't think we should change the semantics of k in the middle of the game.

Given the note and 1 and 2 above, do you agree that code was already good as it was? You would then only have to bring 3D in line.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just for reference, we can easily calculate that in python with scipy (even though I guess we won't include that as a dependency, we can put the code in the comment to regenerate the values):

def pct_to_sigma(pct, dof):
    return np.sqrt(scipy.stats.chi2.ppf(pct / 100., df=dof))

def sigma_to_pct(sigma, dof):
    return scipy.stats.chi2.cdf(sigma**2, df=dof) * 100.

for dim in range(0, 4):
    print("{}D".format(dim), end="")
    for n_sigma in range(1, 6):
        if dim == 0: print("\t    {}    ".format(n_sigma), end="")
        else: print("\t{:.5f}".format(sigma_to_pct(n_sigma, dim)), end="")
    print()
0D	    1    	    2    	    3    	    4    	    5    
1D	68.26895	95.44997	99.73002	99.99367	99.99994
2D	39.34693	86.46647	98.88910	99.96645	99.99963
3D	19.87480	73.85359	97.07091	99.88660	99.99846

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cool, this code is more general than the one I put in plot.py (which I assume you looked at) and it agrees with it in the 2D case. So I propose:

  • undo the changes to 2D plotting
  • we keep k=5 for 2D
  • k == sigma, not variance
  • use k = 5 for 3D as well, which is 99.99846 of all probability.
  • make sure it multiplies sigmas, not variances.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just for reference, we can easily calculate that in python with scipy (even though I guess we won't include that as a dependency, we can put the code in the comment to regenerate the values):

Actually, maybe we should include this code, and comment on some of its output values, just like @senselessDev did below. I think that it makes sense for users to specify an integer number of standard deviations, and we can document in 2D and 3D how much percentage 1,2,3,4,5 stddevs corresponds to. k can be an input argument defaulting to 5. But we could provide another optional argument, percentage that, when given, calculates k using pct_to_sigma.

0D 1 2 3 4 5
1D 68.26895 95.44997 99.73002 99.99367 99.99994
2D 39.34693 86.46647 98.88910 99.96645 99.99963
3D 19.87480 73.85359 97.07091 99.88660 99.99846

U, S, _ = np.linalg.svd(P)

radii = k * np.sqrt(S)
Expand All @@ -113,7 +114,14 @@ def plot_point2_on_axes(axes,
linespec: str,
P: Optional[np.ndarray] = None) -> None:
"""
Plot a 2D point on given axis `axes` with given `linespec`.
Plot a 2D point and its corresponding uncertainty ellipse on given axis
`axes` with given `linespec`.

Based on Stochastic Models, Estimation, and Control Vol 1 by Maybeck,
page 366
For the 2D case:
k = 2.296 corresponds to 1 std, 68.26% of all probability
k = 11.820 corresponds to 3 std, 99.74% of all probability
Copy link
Contributor

@senselessDev senselessDev Jan 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pct_to_sigma(sigma_to_pct(1, 1), 2) ** 2 --> 2.295748928898636
pct_to_sigma(sigma_to_pct(3, 1), 2) ** 2 --> 11.829158081900795

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same


Args:
axes (matplotlib.axes.Axes): Matplotlib axes.
Expand All @@ -125,16 +133,15 @@ def plot_point2_on_axes(axes,
if P is not None:
w, v = np.linalg.eig(P)

# "Sigma" value for drawing the uncertainty ellipse. 5 sigma corresponds
# to a 99.9999% confidence, i.e. assuming the estimation has been
# computed properly, there is a 99.999% chance that the true position
# of the point will lie within the uncertainty ellipse.
k = 5.0
# Scaling value for the uncertainty ellipse, we multiply by 2 because
# matplotlib takes the diameter and not the radius of the major and
# minor axes of the ellipse.
k = 2*np.sqrt(11.820)

angle = np.arctan2(v[1, 0], v[0, 0])
e1 = patches.Ellipse(point,
np.sqrt(w[0] * k),
np.sqrt(w[1] * k),
np.sqrt(w[0]) * k,
np.sqrt(w[1]) * k,
np.rad2deg(angle),
fill=False)
axes.add_patch(e1)
Expand Down Expand Up @@ -178,6 +185,12 @@ def plot_pose2_on_axes(axes,
"""
Plot a 2D pose on given axis `axes` with given `axis_length`.

Based on Stochastic Models, Estimation, and Control Vol 1 by Maybeck,
page 366
For the 2D case:
k = 2.296 corresponds to 1 std, 68.26% of all probability
k = 11.820 corresponds to 3 std, 99.74% of all probability
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same?


Args:
axes (matplotlib.axes.Axes): Matplotlib axes.
pose: The pose to be plotted.
Expand Down Expand Up @@ -205,13 +218,12 @@ def plot_pose2_on_axes(axes,

w, v = np.linalg.eig(gPp)

# k = 2.296
k = 5.0
k = 2*np.sqrt(11.820)

angle = np.arctan2(v[1, 0], v[0, 0])
e1 = patches.Ellipse(origin,
np.sqrt(w[0] * k),
np.sqrt(w[1] * k),
np.sqrt(w[0]) * k,
np.sqrt(w[1]) * k,
np.rad2deg(angle),
fill=False)
axes.add_patch(e1)
Expand Down