椭圆的协方差矩阵

Question

我一直在努力解决一个问题。我很惊讶我没能在网上找到任何有用的东西。

我知道，根据椭圆的协方差矩阵的特征值，可以计算椭圆的长轴和轴长。具体如下：

a1 = 2*sqrt(e1)
a2 = 2*sqrt(e2)

其中a1和a2是长轴和短轴，e1和e2是协方差矩阵的特征值。

我的问题是:给定图像椭圆的边缘点(xi,yi)，如何才能找到该椭圆的2×2协方差矩阵？

Answer 1

仅仅通过纯粹的逆向工程(我不再熟悉这种材料)，我就可以做到这一点：

%// Generate circle
R = 189;
t = linspace(0, 2*pi, 1000).';
x = R*cos(t);
y = R*sin(t);

%// Find the radius?
[~,L] = eig( cov([x,y]) );

%// ...hmm, seems off by a factor of sqrt(2)
2*sqrt( diag(L) )        

%// so it would come out right when I'd include a factor of 1/2 in the sqrt():
2*sqrt( diag(L)/2 )        

那么，让我们检验一下一般椭圆的理论：

%// Random radii
a1 = 1000*rand;
a2 = 1000*rand;

%// Random rotation matrix
R = @(a)[
    +cos(a) +sin(a); 
    -sin(a) +cos(a)];

%// Generate pionts on the ellipse 
t = linspace(0, 2*pi, 1000).';
xy = [a1*cos(t)  a2*sin(t);] * R(rand);

%// Find the deviation from the known radii
%// (taking care of that the ordering may be different)
[~,L] = eig(cov(xy));
min(abs(1-bsxfun(@rdivide, 2*sqrt( diag(L)/2 ), [a1 a2; a2 a1])),[],2)

总是能得到一些小得可以接受的东西。

因此，似乎是可行的:)有人能证实这确实是正确的吗？

Answer 2

为了扩展Rody的答案，固体椭圆的协方差矩阵有lambda_i = r_i^2/4给出的特征值。这引出了OP的r_i = 2*sqrt(lambda_i)方程。

对于(非固态)椭圆，如OP的情况，特征值是固体情况的两倍：lambda_i = r_i^2/2，导致r_i = sqrt(2*lambda_i) (它等于Rody的2*sqrt(lambda_i/2))。

我无法直接找到这方面的参考，但是协方差矩阵的数学与惯性矩的数学是相同的。维基百科 --您可以看到“圆形环”和“实心磁盘”的情况，它们的差异与2的系数相同。

以下是罗迪的测试的改编，既做实心的也做非实心的案例：

% Radius to test with
r = rand(1,2);

% Random rotation matrix
R = @(a)[+cos(a) +sin(a); 
         -sin(a) +cos(a)];

% Generate pionts on the ellipse
N = 1000;
t = linspace(0, 2*pi, N).';
xy = r.*[cos(t),sin(t)] * R(rand);
% Compute radii, compare to known radii
L = eig(cov(xy));
r_ = sqrt(2*L)';
err = max(abs(1 - sort(r_) ./ sort(r)))

% Generate points in the ellipse (SOLID CASE)
N = 10000;
t = 2*pi * rand(N,1);
xy = r .* sqrt(rand(N,1)) .* [cos(t),sin(t)] * R(rand);
% Compute radii, compare to known radii
L = eig(cov(xy));
r_ = 2*sqrt(L)';
err_solid = max(abs(1 - sort(r_) ./ sort(r)))

如果运行这段代码，您将看到1e-3和~6e-3的错误(对于实例，我会生成更多的点，因为这个区域需要更多的点来进行足够密集的采样；点越多，误差就越小)。