Corsario 在评论中提供了一个很好的解决方案：使用核密度函数来测试是否包含在水平集中。

该问题的另一种解释是，它要求一个程序来测试是否包含在由数据的二元正态近似创建的椭圆中。首先，让我们生成一些看起来像问题中的插图的数据：

library(mvtnorm) # References rmvnorm()
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))

椭圆由数据的一阶和二阶矩确定：

center <- apply(p, 2, mean)
sigma <- cov(p)

该公式需要对方差-协方差矩阵求逆：

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

椭圆“高度”函数是二元正态密度的对数的负数：

ellipse <- function(s,t) {u<-c(s,t)-center; u %*% sigma.inv %*% u / 2}

（我忽略了一个等于 $\log(2\pi\sqrt{\det(\Sigma)})$ 的附加常数。）

为了测试这一点，让我们画出它的一些轮廓。这需要在 x 和 y 方向上生成点网格：

n <- 50
x <- (0:(n-1)) * (500000/(n-1))
y <- (0:(n-1)) * (50000/(n-1))

计算此网格的高度函数并绘制它：

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
contour(x,y,matrix(z,n,n), levels=(0:10), col = terrain.colors(11), add=TRUE)

显然它有效。因此，确定点 $(s,t)$ 是否位于 $c$ 水平的椭圆轮廓内的测试是

ellipse(s,t) <= c

Mathematica以同样的方式完成工作：计算数据的方差-协方差矩阵，将其反转，构造ellipse函数，然后一切就绪。

使用R 包的ellipse()功能，该图很简单：mixtools

library(mixtools)
library(mvtnorm) 
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 

第一种方法

您可以在 Mathematica 中尝试这种方法。

让我们生成一些双变量数据：

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

然后我们需要加载这个包：

Needs["MultivariateStatistics`"]

现在：

ellPar=EllipsoidQuantile[data, {0.9}]

给出定义 90% 置信椭圆的输出. 您从此输出中获得的值采用以下格式：

{Ellipsoid[{x1, x2}, {r1, r2}, {{d1, d2}, {d3, d4}}]}

x1 和 x2 指定椭圆居中的点，r1 和 r2 指定半轴半径，d1、d2、d3 和 d4 指定对齐方向。

您还可以绘制此图：

Show[{ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio -> 1],  Graphics[EllipsoidQuantile[data, 0.9]]}]

椭圆的一般参数形式为：

ell[t_, xc_, yc_, a_, b_, angle_] := {xc + a Cos[t] Cos[angle] - b Sin[t] Sin[angle],
    yc + a Cos[t] Sin[angle] + b Sin[t] Cos[angle]}

你可以用这种方式绘制它：

ParametricPlot[
    ell[t, ellPar[[1, 1, 1]], ellPar[[1, 1, 2]], ellPar[[1, 2, 1]], ellPar[[1, 2, 2]],
    ArcTan[ellPar[[1, 3, 1, 2]]/ellPar[[1, 3, 1, 1]]]], {t, 0, 2 \[Pi]},
    PlotRange -> {{0, 100}, {0, 100}}]

您可以根据纯几何信息进行检查：如果椭圆中心 (ellPar[[1,1]]) 与您的数据点之间的欧几里得距离大于椭圆中心与边界之间的距离椭圆（显然，与您的点所在的方向相同），则该数据点位于椭圆之外。

第二种方法

这种方法基于平滑的内核分布。

这些是以与您的数据类似的方式分布的一些数据：

data1 = RandomVariate[BinormalDistribution[{.3, .7}, {.2, .3}, .8], 500];
data2 = RandomVariate[BinormalDistribution[{.6, .3}, {.4, .15}, .8], 500];
data = Partition[Flatten[Join[{data1, data2}]], 2];

我们在这些数据值上获得了平滑的内核分布：

skd = SmoothKernelDistribution[data];

我们为每个数据点获得一个数值结果：

eval = Table[{data[[i]], PDF[skd, data[[i]]]}, {i, Length[data]}];

我们固定一个阈值并选择所有高于此阈值的数据：

threshold = 1.2;
dataIn = Select[eval, #1[[2]] > threshold &][[All, 1]];

在这里，我们得到了区域之外的数据：

dataOut = Complement[data, dataIn];

现在我们可以绘制所有数据：

Show[ContourPlot[Evaluate@PDF[skd, {x, y}], {x, 0, 1}, {y, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, PlotPoints -> 50],
ListPlot[dataIn, PlotStyle -> Darker[Green]],
ListPlot[dataOut, PlotStyle -> Red]]

绿色点是阈值以上的点，红色点是阈值以下的点。

R 包中的ellipse函数ellipse将生成这些椭圆（实际上是一个近似椭圆的多边形）。你可以使用那个椭圆。

实际上可能更容易的是计算您所在点的密度高度，看看它是否高于（椭圆内）或低于（椭圆外）椭圆的轮廓值。函数内部ellipse使用 $\chi^2$ 值来创建椭圆，您可以从那里开始查找要使用的高度。