How do I reorder a regression plot in the desired sequence?
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Akshay Arvind
am 1 Jul. 2024
Kommentiert: Mathieu NOE
am 3 Jul. 2024
I am working on a set of datapoints like so:
x = [270 280 290 300 310 320 330 340 350 0 10 20 30 40 50 60 70 80 90]
y = [10000 9000 5500 2500 900 2500 5500 9000 10000 9000 5500 2500 900 2500 5500 9000 10000]
I have defined a sine^2(x) function with a phase shift of 45 degree, for fitting along these data points, because my data points are shifted that way.
I define x1 = 1:numel(x);
Using set(gca,'xTick',x1,'XTickLabel',x), I know that I can display in the x-axis order 270,...,0,...,90, without which I would get them in the x-axis order 0...90,...,270,...,350.
But, how do I apply that order for the fitting function itself? That doesn't seem to work.
I attach the code for your reference.
function[]=plotdata(filename)
S = load(filename);
C = struct2cell(S);
M = cell2mat(C);
x = M(:,1)
y = M(:,end)
x1 = 1:numel(x);
mean(y)
xlocs = [270 0 90]
%freq = 1/(2*mean(diff(xlocs))) %diff(xlocs)
freq = 1 / (2*mean(diff(xlocs)))
[lb,ub] = bounds(y)
fcn = @(b,x)b(1).*cos(2*pi*x*b(2)+b(3)+(pi/4)).^2+b(4)
B0 = [ub-lb; freq; 0; lb]
myfun = @(b)norm(fcn(b,x) - y);
[B,fv] = fminsearch(myfun,B0)
xv = linspace(max(x),min(x),1000); %A smoother x vector
figure
plot(x, y, '*', 'DisplayName','Data')
hold on
plot(xv, fcn(B,xv), '-r', 'DisplayName','Regression')
hold off
%set(gca,'xTick',x1,'XTickLabel',x)
end
2 Kommentare
Mathieu NOE
am 2 Jul. 2024
why is it important for you to have the data sorted this way ?
that has no influence on the fitting process
Akzeptierte Antwort
Mathieu NOE
am 2 Jul. 2024
hello again
you could do that
NB that the cos^2 is not what I used
A simple sin (or cos if you prefer works better - your y data does not have a squared cos shape
also I added a few missing y points in your data (x and y vectors must be same length)
to make the math a bit simpler for me , I prefer to have x monotonic, so for the 0 to 90 deg segment I added 360 deg
and I simply retrieved again this 360 deg shift to x1 to have the xticklabels as you wanted
hope it helps !
% x = [270 280 290 300 310 320 330 340 350 0 10 20 30 40 50 60 70 80 90];
x = [[270 280 290 300 310 320 330 340 350] [0 10 20 30 40 50 60 70 80 90]+360];
% y = [10000 9000 5500 2500 900 2500 5500 9000 10000 9000 5500 2500 900 2500 5500 9000 10000];
y = [10000 9000 5500 2500 900 2500 5500 9000 10000 9000 5500 2500 900 2500 5500 9000 10000 9000 5500 ];
%%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%
ym = mean(y); % Estimate offset
yu = max(y);
yl = min(y);
yr = (yu-yl)/2; % Estimate Amplitude
% Find zero-crossings
zt = find_zc(x,y,ym);
per = mean(diff(zt)); % Estimate period
fre = 1/per; % Estimate FREQUENCY
phe = -2*pi*zt(1)/per; % Estimate phase
% stationnary sinus fit
fit = @(b,t) b(1)*(sin(2*pi*t*b(2) + b(3))) + b(4); % Objective Function to fit
fcn = @(b) norm(fit(b,x) - y); % Least-Squares cost function
sol = fminsearch(fcn, [yr; fre; phe; ym;]); % Minimise Least-Squares
sol(3) = mod(sol(3),2*pi);
yp = fit(sol,x);
% display sinus fit parameters in command window
amplitude = sol(1);
frequency_Hz = sol(2)
phase_offset_rad = sol(3)
DC_offset = sol(4)
% plot
figure(1),
plot(x, y, 'b',zt, ym*ones(size(zt)), 'dr',x, yp, '-g')
legend('data','zc points','model fit');
xlim([x(1) x(end)]);
x1 = x;
ind = x1>=360;
x1(ind) = x1(ind) - 360;
set(gca,'xTick',x,'XTickLabel',num2str(x1'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
2 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Least Squares finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!