如何找到椭圆内的点？

Question

我试图找到椭圆内的点。它不是一个“普通”椭圆，实际上它基于平均值和标准差，这比计算特征值以找到置信区间要容易得多

函数不是我写的，这里是来源 https://matplotlib.org/devdocs/gallery/statistics/confidence_ellipse.html https://carstenschelp.github.io/2018/09/14/Plot_Confidence_Ellipse_001.html

这是代码：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms

x = np.array([21.5,16.3,13.7,20.0,17.4,10.4,16.9,7.0,13.8,15.2,13.8,8.2,18.0,9.4,13.2,7.2,21.2,30.2,13.5,29.8,18.3,20.2,31.1,21.5,29.8,18.0,13.1,24.1,32.5,15.4,16.1,15.0,25.9,3.0,17.0,23.6,17.6,-11.8,22.2,26.6,17.8,20.6,23.0,28.0,25.3,22.1,22.4,16.3,22.0,12.1])
y = np.array([92.4,98.2,97.6,95.9,96.5,92.1,89.6,89.4,89.2,89.4,90.2,86.7,89.5,89.9,90.2,87.6,104.0,87.3,99.4,85.4,92.8,92.0,87.9,96.2,94.1,95.2,95.6,86.3,87.6,89.5,95.0,97.1,93.0,87.8,98.9,98.2,100.1,45.4,92.1,91.6,94.7,93.9,91.4,91.1,95.7,93.8,96.4,94.1,94.0,89.1])

#function obtained from matplotlib documentation
#https://matplotlib.org/devdocs/gallery/statistics/confidence_ellipse.html

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-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)


#implementation
fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.scatter(x,y,s=5)
ellipse = confidence_ellipse(x,y,ax,n_std=2,edgecolor='red')
plt.show()

后来我尝试找到中心坐标和椭圆内的点，如下所示：

ellipse.get_center()
Out:(0,0)
ellipse.contains_point([21.5,92.4])#first points in x,y arrays
Out:False
ellipse.contains_point([0,0])#get_center() result
Out:False

椭圆图工作正常，但我需要椭圆内的每个点坐标。 我做错了什么？我已经检查过类似的问题，但它们也不起作用。

Answer 1

confidence_ellipse 示例函数仅返回一个用于绘图的对象，并且包含点只会告诉您该点是否在椭圆上。

您可能想要的是这样的：

import math
class distribution():
    def __init__(self,x,y):
        self.cov  = np.cov(x, y)        
        self.mean = np.matrix( [np.mean(x), np.mean(y)])
    def dist(self, x,y):
        tmp = np.matrix([x,y])
        diff = self.mean - tmp
        dist = diff * np.linalg.inv(self.cov) * diff.T
        return math.sqrt(dist)
    def is_inside(self, x,y,nstd=2.0):
    if (self.dist(x,y) < nstd):
        return True
    else:
        return False
        

然后你可以这样做：

d = distribution(x,y)    
d.is_inside(24.1,86.3)

返回 true。

然后，对于所有要点：

points = np.array(list(zip(x, y)))

points_in  = list(filter(lambda p: d.is_inside(p[0],p[1]), points))
points_out =  list(filter(lambda p: not d.is_inside(p[0],p[1]), points))
x_in = [ x[0] for x in points_in] 
y_in = [ x[1] for x in points_in] 

x_out = [ x[0] for x in points_out] 
y_out = [ x[1] for x in points_out] 


fig2, ax2 = plt.subplots(1, 1, figsize=(8, 8))
ax2.scatter(x_in,y_in,s=5, facecolor="green")
ax2.scatter(x_out,y_out, s=5, facecolor="red")
ellipse = confidence_ellipse(x,y,ax2,n_std=2,edgecolor='red') # this presupposes that you still have the confidence_ellipse still defined

plt.show()

你的输出应该是这样的： 红点距离超过 2 个标准差，绿点在内部。

Answer 2

您可以在绘图上绘制所有

x,y

标签。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms

x = np.array([21.5,16.3,13.7,20.0,17.4,10.4,16.9,7.0,13.8,15.2,13.8,8.2,18.0,9.4,13.2,7.2,21.2,30.2,13.5,29.8,18.3,20.2,31.1,21.5,29.8,18.0,13.1,24.1,32.5,15.4,16.1,15.0,25.9,3.0,17.0,23.6,17.6,-11.8,22.2,26.6,17.8,20.6,23.0,28.0,25.3,22.1,22.4,16.3,22.0,12.1])
y = np.array([92.4,98.2,97.6,95.9,96.5,92.1,89.6,89.4,89.2,89.4,90.2,86.7,89.5,89.9,90.2,87.6,104.0,87.3,99.4,85.4,92.8,92.0,87.9,96.2,94.1,95.2,95.6,86.3,87.6,89.5,95.0,97.1,93.0,87.8,98.9,98.2,100.1,45.4,92.1,91.6,94.7,93.9,91.4,91.1,95.7,93.8,96.4,94.1,94.0,89.1])

#function obtained from matplotlib documentation
#https://matplotlib.org/devdocs/gallery/statistics/confidence_ellipse.html

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-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)


#implementation
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
ax.scatter(x,y,s=15)
ellipse = confidence_ellipse(x,y,ax,n_std=2,edgecolor='red')

# zip joins x and y coordinates in pairs
for x,y in zip(x,y):

    label = f"({x},{y})"
    #label = "{:.2f}".format(y) # plot just y-value of the point
    # print(label) # uncomment if you want to print the points for reference 
    plt.annotate(label, # this is the text
                 (x,y), # this is the point to label
                 textcoords="offset points", # how to position the text
                 xytext=(0,10), # distance from text to points (x,y)
                 ha='center') # horizontal alignment can be left, right or center
              
plt.show()

P.S.：您需要相应地调整您的

xytext

。这正是图中绘制的

label

点。

您也可以

print

这些值供您参考。只需输入

print(label)

，它就会为您打印所有分数。

(21.5,92.4)
(16.3,98.2)
(13.7,97.6)
(20.0,95.9)
(17.4,96.5)
(10.4,92.1)
(16.9,89.6)
(7.0,89.4)
(13.8,89.2)
(15.2,89.4)
(13.8,90.2)
(8.2,86.7)
(18.0,89.5)
(9.4,89.9)
(13.2,90.2)
(7.2,87.6)
(21.2,104.0)
(30.2,87.3)
(13.5,99.4)
(29.8,85.4)
(18.3,92.8)
(20.2,92.0)
(31.1,87.9)
(21.5,96.2)
(29.8,94.1)
(18.0,95.2)
(13.1,95.6)
(24.1,86.3)
(32.5,87.6)
(15.4,89.5)
(16.1,95.0)
(15.0,97.1)
(25.9,93.0)
(3.0,87.8)
(17.0,98.9)
(23.6,98.2)
(17.6,100.1)
(-11.8,45.4)
(22.2,92.1)
(26.6,91.6)
(17.8,94.7)
(20.6,93.9)
(23.0,91.4)
(28.0,91.1)
(25.3,95.7)
(22.1,93.8)
(22.4,96.4)
(16.3,94.1)
(22.0,94.0)
(12.1,89.1)

Answer 3

您可以替换原始函数的最后一行 # https://matplotlib.org/stable/gallery/statistics/confidence_ellipse.html#the-plotting-function-itself

...
return ax.add_patch(ellipse)

这两条新线：

...
centroid=[mean_x,mean_y] # mean_x,mean_y are calculated in the function itself and it is the center of the ellipse
return ax.add_patch(ellipse), centroid

然后你像这样调用该函数：

ellipse, centroid = confidence_ellipse(x=x, y=y, ax=ax, edgecolor='red')

它将返回以下内容：

centroid

[18.110000000000003, 91.82200000000003]