Linear data fitting
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi everyone,
So, I have two sets of data, and plotting them against each other I get something like this:
Easily, it is possible to identify 3 groups of data that would be fitted by 3 linear regressions. The problem is that the data points don't follow an order in witch it is possible just to break the arrays in 3 parts and get 3 different data sets. The data points are all mixed up, and there are even points that don't fit in any of these linear fitting groups.
My question is: what is the easiest way for me to break the data set in smaller segments that will be fitted by a linear regression, ignoring the data points that will not fit in any group?
Thank you all in advance.
0 commentaires
Réponse acceptée
Teja Muppirala
le 29 Sep 2011
The problem is that your data is scaled very differently than the random data. If you normalize the data to a better scale first, then fit lines, and then undo the normalization, everything works fine. Convergence to the right answer depends on initial conditions, so I just have it running in an infinite loop until the user stops it. You could certainly implement some check to stop automatically.
thisfig = figure;
load DAT_SET %<----- Run it with your data
plot(x,y,'*');
%Rescale to [0 1]
xmin = min(x); xrng = max(x)-min(x);
ymin = min(y); yrng = max(y)-min(y);
xscl = (x-xmin)/(xrng);
yscl = (y-ymin)/(yrng);
% Sum of squares seems to work faster/better than sum of absolute differences
Pfun = @(P) sum(min(abs(bsxfun(@minus,[P(1)*xscl + P(2) ; P(3)*xscl + P(4) ; P(5)*xscl+ P(6)],yscl)).^2));
options = optimset('display','off','Largescale','off');
options.MaxIter = 1e6;
options.MaxFunEvals = 1e6;
% Runs forever. Close the window when you want it to stop.
fmin = Inf;
hplots = [];
n = 0;
while ishandle(thisfig)
n = n+1;
[P0,f] = fminunc(Pfun,randn(1,6),options);
if f < fmin
fmin = f;
P_opt = P0;
delete(hplots);
xplot = [xmin xmin+xrng];
hold on;
% Convert back to the original coordinates
P_opt([2 4 6]) = -P_opt([1 3 5])*(yrng/xrng)*xmin + P_opt([2 4 6])*(yrng) + ymin;
P_opt([1 3 5]) = P_opt([1 3 5])*(yrng/xrng);
hplots(1) = plot(xplot,P_opt(1)*xplot + P_opt(2),'r');
hplots(2) = plot(xplot,P_opt(3)*xplot + P_opt(4),'r');
hplots(3) = plot(xplot,P_opt(5)*xplot + P_opt(6),'r');
end
title(['Close this figure when you are happy. Iteration number: ' num2str(n)]);
drawnow;
end
% Replot
figure(thisfig);
hold all;
plot(x,y,'*');
plot(xplot,P_opt(1)*xplot + P_opt(2),'r');
plot(xplot,P_opt(3)*xplot + P_opt(4),'r');
plot(xplot,P_opt(5)*xplot + P_opt(6),'r');
Plus de réponses (3)
Fangjun Jiang
le 26 Sep 2011
Certainly this requires some intelligence. I don't see a simple algorithm that can achieve it. What I suggest is to combine programming with your visual inspection.
You've seen the figure so you can roughly find out the fitting for those three straight lines. Use those three lines to calculate the errors for each data point. If the error is large for a particular line, eliminate that data point for that line.
With that, you can separate all the data points into three groups and then you can do the accurate fitting.
0 commentaires
Teja Muppirala
le 27 Sep 2011
One possible way, if you have access to the functionality in the Optimization Toolbox, is to treat this as an optimization problem. Find 6 parameters (2 for each line) that minimize the error amongst the points.
Try it out a few times using this script, it works great!
figure;
x1 = randn(1,100);
x2 = 1+randn(1,100);
x3 = 2+randn(1,100);
Ptrue = 3*randn(1,6);
y1 = Ptrue(1)*x1+Ptrue(2);
y2 = Ptrue(3)*x2+Ptrue(4);
y3 = Ptrue(5)*x3+Ptrue(6);
y1(1:10:end) = randn(size(y1(1:10:end))); % Add noise
plot([x1 x2 x3],[y1 y2 y3],'*')
x = [x1 x2 x3];
y = [y1 y2 y3];
% Use sum of absolute differences (not squares) for robustness
Pfun = @(P) sum(min(abs(bsxfun(@minus,[P(1)*x + P(2) ; P(3)*x + P(4) ; P(5)*x + P(6)],y))));
options = optimset('display','off','Largescale','off');
options.MaxIter = 1e6;
options.MaxFunEvals = 1e6;
% Do it 10 times just to make sure
fmin = Inf;
for n = 1:10
[P0,f] = fminunc(Pfun,randn(1,6),options);
if f < fmin
fmin = f;
P_opt = P0;
end
end
xplot = [min(x), max(x)];
hold on;
plot(xplot,P_opt(1)*xplot + P_opt(2),'r')
plot(xplot,P_opt(3)*xplot + P_opt(4),'r')
plot(xplot,P_opt(5)*xplot + P_opt(6),'r')
Another approach, if you are familiar with image processing, is to use a Hough transform to try and find those lines. For only 3 lines, I think the optimization approach is much easier and more accurate though.
Richard Willey
le 27 Sep 2011
From my perspective, the easiest way to solve this one is
- Fit a linear model to the complete data set
- Apply a clustering algorithm to the residuals
Check out the following
Generate some data
X1 = 1:100;
X1 = X1';
Y1 = 5*X1;
X2 = 1:100;
X2 = X2'+ 1/3;
Y2 = 5* X2 + 50;
X3 = 1:100;
X3 = X3'+ 2/3;
Y3 = 5* X2 + 100;
X = vertcat(X1,X2,X3);
Y = vertcat(Y1,Y2,Y3);
scatter(X,Y, '.')
%%Apply KMeans Clustering to the residuals
B = regress(Y,X);
YHat = X*B;
resid = Y - YHat;
index = kmeans(resid,3);
%%Generate three regression models
for i = 1:3
B(i) = regress(Y(index == i), X(index == i));
end
0 commentaires
Voir également
Catégories
En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!