确定数据的回归系数 - MATLAB

Question

我正在做一个涉及科学计算的项目。以下是我经过一些实验得到的三个变量及其值。

还有一个包含三个未知数的方程，a、b 和 c：

x=(a+0.98)/y+(b+0.7)/z+c

如何使用上述方法获取 a,b,c 的值？这在 MATLAB 中可行吗？

Answer 1

这听起来像是一个回归问题。假设测量中无法解释的误差是高斯分布的，你可以通过least squares找到参数。基本上，您必须重写方程，以便将其变为ma + nb + oc = p 的形式，然后您有 6 个具有 3 个未知数的方程 (a, b, c)，这些参数可以通过最小二乘法优化找到。因此，通过一些代数，我们得到：

za + yb + yzc = xyz - 0.98z - 0.7z

因此，m = z, n = y, o = yz, p = xyz - 0.98z - 0.7z。我将把它留给你作为练习，以验证我的代数是否正确。然后可以形成矩阵方程：

Ax = d

我们将有 6 个方程，我们想求解 x 其中x = [a b c]^{T}。要求解x，您可以使用所谓的pseudoinverse 来检索参数，以最大程度地减少真实输出与使用相同输入数据时由这些参数生成的输出之间的误差。

换句话说：

x = A^{+}d

A^{+} 是矩阵A 的伪逆矩阵，是矩阵向量乘以向量d。

为了将我们的想法转化为代码，我们将定义输入数据，形成A 矩阵和d 向量，其中它们之间共享的每一行都是一个方程，然后使用伪逆来找到我们的参数。您可以使用ldivide (\) 运算符来完成这项工作：

%// Define x y and z
x = [9.98; 8.3; 8.0; 7; 1; 12.87];
y = [7.9; 7.5; 7.4; 6.09; 0.9; 11.23];
z = [7.1; 5.6; 5.9; 5.8; -1.8; 10.8];

%// Define A matrix
A = [z y y.*z];
%// Define d vector
d = x.*y.*z - 0.98*z - 0.7*z;

%// Find parameters via least-squares
params = A\d;

params存储参数a、b和c，我们得到：

params =

  -37.7383
  -37.4008
   19.5625

---

如果您想仔细检查这些值的接近程度，您只需在帖子中使用上述表达式并与x 中的每个值进行比较：

a = params(1); b = params(2); c = params(3);
out = (a+0.98)./y+(b+0.7)./z+c;
disp([x out])

9.9800    9.7404
8.3000    8.1077
8.0000    8.3747
7.0000    7.1989
1.0000   -0.8908
12.8700   12.8910

您可以看到它并不完全接近，但您获得的参数在 最小二乘 错误意义上是最好的。

---

奖金 - 与 RANSAC 配合

您可以看到一些预测值（输出中的右列）比其他预测值更多。这是因为我们使用了您数据中的所有点来找到合适的模型。一种用于最小化错误并提高模型估计稳健性的技术是使用称为 RANSAC 或 RANdom SAmple C强烈>共识。 RANSAC 背后的基本方法是，对于一定数量的迭代，您获取数据并随机抽取找到模型所需的最少点。找到此模型后，如果要使用这些参数来描述数据，就会发现总体错误。你不断地随机选择点，找到你的模型，找到错误，产生最少错误的迭代将是你定义整个模型的参数。

正如您在上面看到的，我们可以定义的一个错误是真实x 点和预测x 点之间的绝对差之和。还有许多其他度量，例如平方误差之和，但现在让我们坚持一些简单的方法。如果你看一下上面的公式，我们至少需要 三个方程 来定义 a、b 和 c，因此对于每次迭代，我们会随机选择三个点无需替换我可能会添加，找到我们的模型，确定错误，然后继续迭代并找到错误最少的参数。

因此，您可以像这样编写 RANSAC 算法：

%// Define cost and number of iterations
cost = Inf;
iterations = 50;

%// Set seed for reproducibility
rng(123);

%// Define x y and z
x = [9.98; 8.3; 8.0; 7; 1; 12.87];
y = [7.9; 7.5; 7.4; 6.09; 0.9; 11.23];
z = [7.1; 5.6; 5.9; 5.8; -1.8; 10.8];

for idx = 1 : iterations
    %// Determine where we would need to sample
    ind = randperm(numel(x), 3);

    xs = x(ind); ys = y(ind); zs = z(ind); %// Sample

    %// Define A matrix
    A = [zs ys ys.*zs];
    %// Define d vector
    d = xs.*ys.*zs - 0.98*zs - 0.7*zs;

    %// Find parameters via least-squares
    params = A\d;

    %// Determine error
    a = params(1); b = params(2); c = params(3);
    out = (a+0.98)./y+(b+0.7)./z+c;
    err = sum(abs(x - out));

    %// If error produced is less than current error
    %// then save parameters
    if err < cost
        cost = err;
        final_params = params;
    end
end

---

当我运行上面的代码时，我得到了我的参数：

final_params =

  -38.1519
  -39.1988
   19.7472

将此与我们的x 进行比较，我们得到：

a = final_params(1); b = final_params(2); c = final_params(3);
out = (a+0.98)./y+(b+0.7)./z+c;
disp([x out])

9.9800    9.6196
8.3000    7.9162
8.0000    8.1988
7.0000    7.0057
1.0000   -0.1667
12.8700   12.8725

如您所见，数值有所提高——尤其是第四点和第六点……并将其与之前的版本进行比较：

9.9800    9.7404
8.3000    8.1077
8.0000    8.3747
7.0000    7.1989
1.0000   -0.8908
12.8700   12.8910

您可以看到第二个值比以前的版本更差，但其他数字更接近真实值。

---

玩得开心！