如何从任一侧延伸 2 PolyFit 的线相交并获得组合拟合线

How to extend the line of 2 PolyFit from either side to intersect and get a combined fit line

我正在尝试从任一侧(应该相交)获得由两条线性 polyfit 制成的组合拟合线,这是拟合线的图片:

我正在尝试使两条拟合(蓝色)线相交并生成一条组合拟合线,如下图所示:

请注意,波峰可以出现在任何地方,所以我不能假设它位于中心。

这是创建第一个图的代码:

xdatPart1 = R;
zdatPart1 = z;

n = 3000;
ln = length(R);

[sX,In] = sort(R,1);

sZ = z(In);       

xdatP1 = sX(1:n,1);
zdatP1 = sZ(1:n,1);

n2 = ln - 3000;

xdatP2 = sX(n2:ln,1);
zdatP2 = sZ(n2:ln,1);

pp1 = polyfit(xdatP1,zdatP1,1);
pp2 = polyfit(xdatP2,zdatP2,1);

ff1 = polyval(pp1,xdatP1);
ff2 = polyval(pp2,xdatP2);

xDat = [xdatPart1];
zDat = [zdatPart1];

axes(handles.axes2);
cla(handles.axes2);
plot(xdatPart1,zdatPart1,'.r');
hold on
plot(xdatP1,ff1,'.b');
plot(xdatP2,ff2,'.b');
xlabel(['R ',units]);
ylabel(['Z ', units]);
grid on
hold off

下面是没有曲线拟合工具箱的粗略实现。虽然代码应该是不言自明的,但这里是算法的大纲:

  1. 我们生成了一些数据。
  2. 我们通过平滑数据并找到最大值的位置来估计交点。
  3. 我们在估计的交点的每一侧拟合一条线。
  4. 我们使用拟合方程计算拟合线的交点。
  5. 我们使用 mkpp 构造 "evaluateable" 分段多项式的函数句柄。
  6. 输出 ppfunc 是 1 个变量的函数句柄,您可以像 any regular function 一样使用它。

现在,这个解决方案在任何意义上都不是最优(例如 MMSE、LSQ 等),但正如您将在与 MATLAB 工具箱的结果比较中看到的那样,还不错!

function ppfunc = q40160257
%% Define the ground truth:
center_x = 6 + randn(1);
center_y = 78.15 + 0.01 * randn(1);
% Define a couple of points for the left section
leftmost_x = 0;
leftmost_y = 78.015 + 0.01 * randn(1);
% Define a couple of points for the right section
rightmost_x = 14.8;
rightmost_y = 78.02 + 0.01 * randn(1);
% Find the line equations:
m1 = (center_y-leftmost_y)/(center_x-leftmost_x);
n1 = getN(leftmost_x,leftmost_y,m1);
m2 = (rightmost_y-center_y)/(rightmost_x-center_x);
n2 = getN(rightmost_x,rightmost_y,m2);
% Print the ground truth:
fprintf(1,'The line equations are: {y1=%f*x+%f} , {y2=%f*x+%f}\n',m1,n1,m2,n2)
%% Generate some data:
NOISE_MAGNITUDE = 0.002;
N_POINTS_PER_SIDE = 1000;
x1 = linspace(leftmost_x,center_x,N_POINTS_PER_SIDE);
y1 = m1*x1+n1+NOISE_MAGNITUDE*randn(1,numel(x1));
x2 = linspace(center_x,rightmost_x,N_POINTS_PER_SIDE);
y2 = m2*x2+n2+NOISE_MAGNITUDE*randn(1,numel(x2));
X = [x1 x2(2:end)]; Y = [y1 y2(2:end)];
%% See what we have:
figure(); plot(X,Y,'.r'); hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimating the intersection point:
MOVING_AVERAGE_PERIOD = 10; % Play around with this value.
smoothed_data = conv(Y, ones(1,MOVING_AVERAGE_PERIOD)/MOVING_AVERAGE_PERIOD, 'same');
plot(X, smoothed_data, '-b'); ylim([floor(leftmost_y*10) ceil(center_y*10)]/10);
[~,centerInd] = max(smoothed_data);
fprintf(1,'The real intersection is at index %d, the estimated is at %d.\n',...
           N_POINTS_PER_SIDE, centerInd);
%% Fitting a polynomial to each side:
p1 = polyfit(X(1:centerInd),Y(1:centerInd),1);
p2 = polyfit(X(centerInd+1:end),Y(centerInd+1:end),1);
[x_int,y_int] = getLineIntersection(p1,p2);
plot(x_int,y_int,'sg');

pp = mkpp([X(1) x_int X(end)],[p1; (p2 + [0 x_int*p2(1)])]);
ppfunc = @(x)ppval(pp,x);
plot(X, ppfunc(X),'-k','LineWidth',3)
legend('Original data', 'Smoothed data', 'Computed intersection',...
       'Final piecewise-linear fit');
grid on; grid minor;     

%% Comparison with the curve-fitting toolbox:
if license('test','Curve_Fitting_Toolbox')
  ft = fittype( '(x<=-(n2-n1)/(m2-m1))*(m1*x+n1)+(x>-(n2-n1)/(m2-m1))*(m2*x+n2)',...
                'independent', 'x', 'dependent', 'y' );
  opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
  % Parameter order: m1, m2, n1, n2:
  opts.StartPoint = [0.02 -0.02 78 78];
  fitresult = fit( X(:), Y(:), ft, opts);
  % Comparison with what we did above:
  fprintf(1,[...
    'Our solution:\n'...
    '\tm1 = %-12f\n\tm2 = %-12f\n\tn1 = %-12f\n\tn2 = %-12f\n'...
    'Curve Fitting Toolbox'' solution:\n'...
    '\tm1 = %-12f\n\tm2 = %-12f\n\tn1 = %-12f\n\tn2 = %-12f\n'],...
    m1,m2,n1,n2,fitresult.m1,fitresult.m2,fitresult.n1,fitresult.n2);    
end

%% Helper functions:
function n = getN(x0,y0,m)
% y = m*x+n => n = y0-m*x0;
n = y0-m*x0;

function [x_int,y_int] = getLineIntersection(p1,p2)
% m1*x+n1 = m2*x+n2 => x = -(n2-n1)/(m2-m1)
x_int = -(p2(2)-p1(2))/(p2(1)-p1(1));
y_int = p1(1)*x_int+p1(2);

结果(样本运行):

Our solution:
    m1 = 0.022982    
    m2 = -0.011863   
    n1 = 78.012992   
    n2 = 78.208973   
Curve Fitting Toolbox' solution:
    m1 = 0.022974    
    m2 = -0.011882   
    n1 = 78.013022   
    n2 = 78.209127   

放大交叉路口: