Comment étendre la ligne de 2 PolyFit de chaque côté à l'intersection et obtenir une ligne d'ajustement combiné

Question

J'essaie d'obtenir une ligne d'ajustement combinée faite de deux polyfit linéaire de chaque côté (devrait se croiser), voici l'image des lignes de forme :

Je suis en train de faire les deux lignes en forme (bleu) se croisent et produisent une ligne en forme combinée comme le montre l'image ci-dessous:

Notez que la crête peut arriver n'importe où donc je ne peux pas supposer être au centre.

Voici le code qui crée la première parcelle:

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 

Answer 1

est une implémentation grossière ci-dessous sans boîte à outils d'ajustement de courbe. Bien que le code devrait être explicite, voici un aperçu de l'algorithme:

1. Nous générons des données.
2. Nous estimons le point d'intersection en lissant les données et en trouvant l'emplacement de la valeur maximale.
3. Nous ajustons une ligne de chaque côté du point d'intersection estimé.
4. Nous calculons l'intersection des lignes ajustées en utilisant les équations ajustées.
5. Nous utilisons mkpp pour construire une poignée de fonction à un polynôme par morceaux «évaluable».
6. La sortie, ppfunc, est un handle de fonction de 1 variable, que vous pouvez utiliser comme any regular function.

Maintenant, cette solution n'est pas optimale dans un sens (comme MMSE, LSQ, etc.), mais comme vous le verrez dans la comparaison avec le résultat de la boîte à outils de Matlab, ce n'est pas si mal que ça!

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); 

Le résultat (échantillon analysé):

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 

zoomant autour de l'intersection: