Нахождение нечетной точки в наборе данных без использования циклов

Question

Мне дан набор точек (p1,q1) (p2,q2) ... (p20,q20), которые удовлетворяют функции q = 1/(ap + b)^2, за исключением того, что один из них не удовлетворяет данному соотношению .Значения a и b мне не даны.Все, что у меня есть, это два входа p и q в виде массивов.Мне нужно найти индекс точки, которая не удовлетворяет данному соотношению.

Я решил найти способ найти значения a и b, используя первые две пары (p1,q1) и (p2,q2), и проверить, удовлетворяют ли оставшиеся точки функции для найденных значений.a и 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

Давайте выполним некоторые математические операции, чтобы избежать необходимости векторизации или . Самое большее, это оставляет нам полдюжины оценок функций, и нам нужно всего 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 )

Answer 2

Я не очень понимаю, как вы пытались решить эту проблему, и что общего с этим связано syms (т.е. символические переменные), поэтому я покажу вам, как I решит этопроблема.

Поскольку мы, по сути, ищем выброс, мы могли бы также преобразовать проблему во что-то, с чем легче работать.По этой причине вместо использования q как есть я собираюсь инвертировать его: таким образом, мы будем иметь дело с уравнением параболы - что легко.

Далее, знаячто наши точки должны лежать на параболе, мы можем подобрать уравнение параболы (или эквивалентно - найти коэффициенты полинома, описывающего отношение входа к выходу).Полином равен a^2*x^2+2*a*b*x+b^2, поэтому коэффициенты равны {a^2, 2*a*b, b^2}.

Поскольку большинство точек (19 из 20) лежат на одной и той же параболе, выброс всегда будет иметь большую ошибку, которая вне зависимости от того, насколько близко она находится к параболе (в пределах точности машинной точности) - вы можете увидеть крайний пример этого в приведенном ниже коде.

Подгонкапараболы выполняется с использованием полиномиальной интерполяции (см. также: матрица Вандермонда ).

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