在 MATLAB 中将平面拟合到 N 维点

Question

我有一组 N 个 k 维度的点作为大小为 N X k 的矩阵。

如何通过这些点找到最佳拟合线？该线将是 k 维度的平面（超平面）。它将有 k 个系数和一个偏置项。

像 fit 这样的现有函数似乎只能用于 2 维或 3 维的点。

Answer 1

您可以使用主成分分析将超平面（或任何低维仿射 space）拟合到一组 D 维数据。下面是一个将平面拟合到一组 3D 数据的示例。这在 MATLAB documentation 中有更详细的解释，但我试图构建我能构建的最简单的示例。

% generate some random correlated data
D = 3;
mu = zeros(1,D);
sqrt_sig = randn(D);
sigma = sqrt_sig'*sqrt_sig;
% generate 50 points in a D x 50 matrix
X = mvnrnd(mu, sigma, 50)';

% perform PCA
coeff = pca(X');

% The last principal component is normal to the best fit plane and plane goes through mean of X
a = coeff(:,D);
b = -mean(X,2)'*a;

% plane defined by a'*x + b = 0
dist = abs(a'*X+b) / norm(a);
mse = mean(dist.^2)

编辑： 添加了 D = 3 的结果示例图。我在这里利用了其他主成分的正交性。如果你想要它只是为了证明飞机确实非常适合数据，请忽略代码。

% plot in 3D
X0 = bsxfun(@minus,X,mean(X,2));
b1 = coeff(:,1);    b2 = coeff(:,2);
y1 = b1'*X0;        y2 = b2'*X0;
y1_min = min(y1);   y1_max = max(y1);
y1_span = y1_max - y1_min;
y2_min = min(y2);   y2_max = max(y2);
y2_span = y2_max - y2_min;
pad = 0.2;
y1_min = y1_min - pad*y1_span;
y1_max = y1_max + pad*y1_span;
y2_min = y2_min - pad*y2_span;
y2_max = y2_max + pad*y2_span;
[y1_m,y2_m] = meshgrid(linspace(y1_min,y1_max,5), linspace(y2_min,y2_max,5));
grid = bsxfun(@plus, bsxfun(@times,y1_m(:)',b1) + bsxfun(@times,y2_m(:)',b2), mean(X,2));
x = reshape(grid(1,:),size(y1_m));
y = reshape(grid(2,:),size(y1_m));
z = reshape(grid(3,:),size(y1_m));

figure(1); clf(1);
surf(x,y,z,'FaceColor','black','FaceAlpha',0.3,'EdgeAlpha',0.6);
hold on;
plot3(X(1,:),X(2,:),X(3,:),' .');
axis equal;
axis vis3d;

Edit2：当我说 "principal component" 时，我的措辞有点草率（或者完全错误）。我实际上指的是表示主成分的正交基向量。

Answer 2

这是一个更简单的解决方案，只使用 MATLAB's \ operator。我们从定义 k 维度的平面开始：

%  0 = a + x(1) * b(1) + x(2) * b(2) + ... + x(k) * 1
k = 8;
a = randn(1);
b = randn(k-1,1);

（注意我们假设b(k)=1，你总是可以在不改变平面的情况下将平面参数乘以任何值）。

接下来我们在这个平面内生成N个随机点：

N = 1000;
x = rand(N,k-1);
x(:,k) = -(a + x * b);

...抱歉，这不是在平面上生成随机点的最佳方式，但对于这里的演示来说已经足够了。为点添加噪声：

x = x + 0.05*randn(size(x));

为了找到飞机的参数，我们求解方程组

%    a + x(1:k-1) * b == -x(k)

在最小二乘意义上。 a 和 b 是那里的未知数。我们可以将左侧重写为 [1,x(1:k-1)] * [a;b]。如果我们有一个矩阵方程 M*p=v 我们可以通过写 p=M\v:

来求解 p

p = [ones(N,1),x(:,1:k-1)]\(-x(:,k));

disp(['ground truth: [a,b,1] = ',mat2str([a,b',1],3)]);
disp(['estimated   : [a,b,1] = ',mat2str([p',1],3)]);

输出为：

ground truth: [a,b,1] = [-1.35 -1.44 -1.48 1.17 0.226 -0.214 0.234 -1.59 1]
estimated   : [a,b,1] = [-1.41 -1.38 -1.43 1.14 0.219 -0.195 0.221 -1.54 1]

数据集中的噪声越少或点越多，误差当然会越小！

在 MATLAB 中将平面拟合到 N 维点

Fit plane to N dimensional points in MATLAB

matlab

curve-fitting

data-fitting