如何从任一侧延伸2 PolyFit的线以相交并获得组合的拟合线 [英] How to extend the line of 2 PolyFit from either side to intersect and get a combined fit line

问题描述

我正在尝试从两侧的两个线性polyfit(应该相交)获得组合的拟合线，这是拟合线的图片:

I am trying to get combined fit line made from two linear polyfit from either side (should intersect), here is the picture of fit lines:

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

I am trying to make the two fit (blue) lines intersect and produce a combined fit line as shown in the picture below:

请注意，波峰可以在任何地方发生，所以我不能假设它在中心.

Note that the crest can happen anywhere so I cannot assume to be in the center.

以下是创建第一个绘图的代码:

Here is the code that creates the first plot:

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

推荐答案

下面是一个粗略的实现，没有曲线拟合工具箱.尽管代码应该是不言自明的，但这是该算法的概述:

Below's a rough implementation with no curve fitting toolbox. Although the code should be self-explanatory, here's an outline of the algorithm:

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

1. We generate some data.
2. We estimate the intersection point by smoothing the data and finding the location of the maximum value.
3. We fit a line to each side of the estimated intersection point.
4. We compute the intersection of the fitted lines using the fitted equations.
5. We use mkpp to construct a function handle to an "evaluateable" piecewise polynomial.
6. The output, ppfunc, is a function handle of 1 variable, that you can use just like any regular function.

现在，该解决方案在任何意义上都不是最优的(例如MMSE，LSQ等)，但是正如您将与MATLAB工具箱的结果进行比较所看到的那样，还不错！

Now, this solution is not optimal in any sense (such as MMSE, LSQ, etc.) but as you will see in the comparison with the result from MATLAB's toolbox, it's not that bad!

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   

在交叉点处放大:

这篇关于如何从任一侧延伸2 PolyFit的线以相交并获得组合的拟合线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！