在不使用循环的情况下查找数据集中的奇点

Question

我得到一组点 (p1,q1) (p2,q2) ... (p20,q20)，它们满足函数 q = 1/(ap + b)^2，除了 其中一个不满足给定的关系。 a 和 b 的值没有给我。我所拥有的只是两个输入 p 和 q 作为数组。我需要找到不满足给定关系的点的索引。

我继续求解的方法是使用前两对 (p1,q1) 和 (p2,q2) 找到 a 和 b 的值，并检查剩余点是否满足 @987654333 求解值的函数@ 和 b。结果将存储在逻辑矩阵中。 我希望利用逻辑矩阵来挑选奇数对，但无法继续进行。

具体来说，挑战在于利用 MATLAB 中的矢量化来找到奇点，而不是诉诸 for 循环。我认为我必须首先在任何行中搜索唯一的逻辑零。在那种情况下，那个零的列索引会给我带来奇数点。但是，如果所有 4 行中有多个零，则奇数点是前两对中的任何一个。我需要帮助将其转换为 MATLAB 中的高效代码。

请注意，向量p 和q 在下面的代码中被命名为x 和y。

function [res, sol] = findThePair(x, y)

N = length(x);

syms a b
vars = [a,b];
eqns = [y(1) - 1/(a*x(1) + b)^2 == 0; y(2) - 1/(a*x(2) + b)^2];
[solA, solB] = solve(eqns,vars);
sol = [double(solA) double(solB)];    %solution of a & b (total 4 possibilites)

xTest = x(3:end);   % performing check on remaining points
yTest = y(3:end);
res = zeros(4, N-2);    % logical matrix to store the results of equality check

for i = 1:4
    A = sol(i,1); B = sol(i, 2);
    res(i, :) = [yTest == 1./(A*xTest + B).^2]; % perform equality check on remaining points
end

Answer 1

我不太明白你是如何解决这个问题的，syms（即符号变量）与这个有什么关系，所以我将向你展示 我 将如何解决这个问题。

由于我们本质上是在寻找异常值，因此我们不妨将问题转换为更易于处理的问题。出于这个原因，我不按原样使用q，而是将其反转：这样，我们将处理抛物线方程——这很容易。

接下来，知道我们的点应该位于抛物线上，我们可以拟合抛物线方程（或等效地 - 找到描述输入与输出关系的多项式的系数）。多项式为a^2*x^2+2*a*b*x+b^2，因此系数为{a^2, 2*a*b, b^2}。

由于大多数点（20 个点中的 19 个）位于同一抛物线上，因此异常值总是会有较大的误差，这会使其脱颖而出，不管它有多接近到抛物线（在机器精度的限制内）——你可以在下面的代码中看到一个极端的例子。

使用polynomial interpolation 执行抛物线拟合（另请参阅：Vandermonde matrix）。

function I = q55241683()
%% Generate the ground truth:
TRUE_A = 2.3;
TRUE_B = -pi;
IDX_BAD = 5;

p = 1:0.04:1.76;
q = (TRUE_A * p + TRUE_B).^-2;
q(IDX_BAD) = (1-1E-10)*q(IDX_BAD); % notice just how close this is to being valid

%% Visualize dataset:
% figure(); plot(p,q.^-1);

%% Solve
I = findThePair(p, q.^-1);

%% Test
if IDX_BAD == I
  disp('Great success!');
else
  disp('Complete failure!');
end

end

function I = findThePair(x,y)
% Fit a parabola to {x vs. y^-1}
P = x(:).^(2:-1:0)\y(:); %alternatively: P = polyfit(x,y.^-1,2)
% Estimate {a,b} (or {-a,-b})
est_A = sqrt(P(1));
est_B = P(2)/(2*est_A);
% Compute the distances of the points from the fit (residuals), find the biggest:
[~,I] = max( abs(y - (est_A*x + est_B).^2) );
end

Answer 2

让我们先做一些数学运算，以避免需要循环或向量化。最多只剩下六次函数评估，我们只需要 5 分。

q = 1 / (a*p + b)^2
% ->
sqrt(q) * ( a*p + b ) = 1
% ->
a = ( 1 - b*sqrt(q) ) / ( p * sqrt(q) )

% Sub in some points (1 and 2) ->
a1 = ( 1 - b*sqrt(q1) ) / ( p1 * sqrt(q1) )    
a2 = ( 1 - b*sqrt(q2) ) / ( p2 * sqrt(q2) )
% a1 and a2 should be the same ->
( 1 - b*sqrt(q1) ) * ( p2 * sqrt(q2) ) = ( 1 - b*sqrt(q2) ) * ( p1 * sqrt(q1) )
% Rearrange ->
b = ( p2*sqrt(q2) - p1*sqrt(q1) ) / ( (p2-p1)*sqrt(q1)*sqrt(q2) )

我们有两个未知数，a 和 b。我们只需要两点来创建联立方程。我将使用以下逻辑

1. 选择(pm, qm) 和(pn, qn) 与任何m ~= n。

2. 使用上述公式计算a 和b。

3. 测试(pr, qr) 是否与计算出的a 和b 相符。

   • 如果合适，我们知道这三个都必须在曲线上，我们有 a 和 b。

   • 如果不合适，我们知道点 m、n 或 r 是异常值。回到步骤（1），另外两个点，计算的a 和b 必须是正确的，因为我们没有拟合异常值。

这里有一些代码来实现这个：

% Random coeffs, keep things unknown
a = rand*10;
b = rand*10;
% Set up our data
p = 1:20;
q = 1 ./ (a*p + b).^2;
% Create an outlier
q( 3 ) = q( 3 ) + 1;

% Steps as described 

% 1.
p1 = p(1); p2 = p(2);
q1 = q(1); q2 = q(2);

% 2.
bGuess = ( p2*sqrt(q2) - p1*sqrt(q1) ) / ( (p2-p1)*sqrt(q1)*sqrt(q2) );
aGuess = ( 1 - bGuess*sqrt(q1) ) / ( p1 * sqrt(q1) );

% 3.
p3 = p(3);
q3Guess = 1 / ( aGuess*p3 + bGuess )^2;

tol = 1e-7; % Use tolerance rather than == comparison to avoid float issues

if abs( q3Guess - q(3) ) < tol
    % success
    aFit = aGuess;
    bFit = bGuess;
else
    % p1, p2 or p3 is an outlier! Repeat using other points
    % If there's known to be only one outlier, this should give the result
    p1 = p(4); p2 = p(5);
    q1 = q(4); q2 = q(5);
    bFit = ( p2*sqrt(q2) - p1*sqrt(q1) ) / ( (p2-p1)*sqrt(q1)*sqrt(q2) );
    aFit = ( 1 - bFit*sqrt(q1) ) / ( p1 * sqrt(q1) );    
end

% Validate
fprintf( 'a is valid: %d, b is valid: %d\n', abs(a-aFit)<tol, abs(b-bFit)<tol )