How to fit curve using for loop?

Question

AS 2021-10-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1466364-how-to-fit-curve-using-for-loop

评论： Mathieu NOE 2021-10-25

I have a 15 number of signals with data point 64 in each. I want to use loop to fit all the signals with fittype 'gauss2' and plot all the curves with thier fitting. I have written like this but it is showing eroor. Kindly suggest me to resolve this. Thank you.

figure;

for i1=1:64

for j1=1:15

recon1_f(j1)=fit(t(i1),recon_amp2_1(:,j1),'gauss2');

h2{i}=plot(recon1_f(j1),t,recon_amp2_1(:,j1));

ylim([0 0.05]);

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Mathieu NOE 2021-10-4

hi

maybe you should share some data so we can test the code

tx

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2021-10-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1466364-how-to-fit-curve-using-for-loop#answer_801254

在 MATLAB Online 中打开

hello again

in the mean time I created this example (only 5 loops) on dummy data

you can easily expand and adapt it to your own needs

NB you only need one for loop and not two as in your code

hope it helps

clearvars
clc
% dummy data
x1 = [0:1:20];
y1 = [0,0.004,0.008,0.024,0.054,0.112,0.33,0.508,0.712,0.926,1,0.874,0.602,0.404,0.252,0.146,0.074,0.036,0.018,0.004,0];
for ci = 1:5
    % modify x and y range (dummy data generation)
    x = x1*ci;
    y = y1*ci^2 + 0.1*rand(size(y1));
    
    % curve fit using fminsearch
    f = @(a,b,c,x) a.*exp(-(x-b).^2 / c.^2);
    obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
    sol = fminsearch(obj_fun, [max(y),max(x)/2,max(x)/6]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    xx = linspace(min(x),max(x),300);
    y_fit = f(a_sol, b_sol,c_sol, xx);
    yy = interp1(x,y, xx);
    Rsquared = my_Rsquared_coeff(yy,y_fit); % correlation coefficient
    figure(ci)
    plot(xx, y_fit, '-',x,y, 'r .', 'MarkerSize', 40)
    title(['Gaussian Fit / R² = ' num2str(Rsquared) ], 'FontSize', 15)
    ylabel('Intensity (arb. unit)', 'FontSize', 14)
    xlabel('x(nm)', 'FontSize', 14)
    eqn = " y = "+a_sol+ " * exp(-(x - " +b_sol+")² / (" +c_sol+ ")²";
    legend(eqn)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

20 个评论
显示 18更早的评论隐藏 18更早的评论

Mathieu NOE 2021-10-5

在 MATLAB Online 中打开

hello

I adpted the code to your data so it will fit two gaussian curves per data row

try this :

clearvars
clc
% load data
x = csvread('time.csv');
y = csvread('data_che.csv');
for ci = 1:5
    
    % curve fit using fminsearch
    % select first peak data
    ind1 = find(x<1700);
    [x_fit1,y_fit1,sol1] = myfit(x(1,ind1),y(ci,ind1));
    a_sol1 = sol1(1);
    b_sol1 = sol1(2);
    c_sol1 = sol1(3);
    y_int1 = interp1(x(1,ind1),y(ci,ind1), x_fit1);
    
    % select second peak data
    ind2 = find(x>=1700);
    [x_fit2,y_fit2,sol2] = myfit(x(1,ind2),y(ci,ind2));
    a_sol2 = sol2(1);
    b_sol2 = sol2(2);
    c_sol2 = sol2(3);
    y_int2 = interp1(x(1,ind2),y(ci,ind2), x_fit2);
    
    
    Rsquared1 = my_Rsquared_coeff(y_int1,y_fit1); % correlation coefficient
    Rsquared2 = my_Rsquared_coeff(y_int2,y_fit2); % correlation coefficient
    Rsquared_global = my_Rsquared_coeff([y_int1 y_int2],[y_fit1 y_fit2]); % correlation coefficient
    figure(ci)
    plot(x_fit1, y_fit1, '-',x_fit2, y_fit2, '-',x,y(ci,:), 'r .', 'MarkerSize', 25)
%     title(['Gaussian Fit / R1² = ' num2str(Rsquared1) ' / R2² = ' num2str(Rsquared2)], 'FontSize', 15)
    title(['Gaussian Fit / R² = ' num2str(Rsquared_global)], 'FontSize', 15)
    ylabel('Y', 'FontSize', 14)
    xlabel('X', 'FontSize', 14)
    eqn1 = " y = "+a_sol1+ " * exp(-(x - " +b_sol1+")² / (" +c_sol1+ ")²";
    eqn2 = " y = "+a_sol2+ " * exp(-(x - " +b_sol2+")² / (" +c_sol2+ ")²";
    legend(eqn1,eqn2,'data')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function     [x_fit,y_fit,sol] = myfit(x,y)
    
    f = @(a,b,c,x) a.*exp(-(x-b).^2 / c.^2);
    obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
    sol = fminsearch(obj_fun, [max(y),max(x)/2,max(x)/6]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    x_fit = linspace(min(x),max(x),100);
    y_fit = f(a_sol, b_sol,c_sol, x_fit);
    
    
end
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

Mathieu NOE 2021-10-5

在 MATLAB Online 中打开

see below. Now the a,b,c solutions are stored as array like :

    a_sol1(ci) = sol1(1);
    b_sol1(ci) = sol1(2);
    c_sol1(ci) = sol1(3);

full code :

    clearvars
clc
% load data
x = csvread('time.csv');
y = csvread('data_che.csv');
for ci = 1:5
    
    % curve fit using fminsearch
    % select first peak data
    ind1 = find(x<1700);
    [x_fit1,y_fit1,sol1] = myfit(x(1,ind1),y(ci,ind1));
    a_sol1(ci) = sol1(1);
    b_sol1(ci) = sol1(2);
    c_sol1(ci) = sol1(3);
    y_int1 = interp1(x(1,ind1),y(ci,ind1), x_fit1);
    
    % select second peak data
    ind2 = find(x>=1700);
    [x_fit2,y_fit2,sol2] = myfit(x(1,ind2),y(ci,ind2));
    a_sol2(ci) = sol2(1);
    b_sol2(ci) = sol2(2);
    c_sol2(ci) = sol2(3);
    y_int2 = interp1(x(1,ind2),y(ci,ind2), x_fit2);
    
    
    Rsquared1 = my_Rsquared_coeff(y_int1,y_fit1); % correlation coefficient
    Rsquared2 = my_Rsquared_coeff(y_int2,y_fit2); % correlation coefficient
    Rsquared_global = my_Rsquared_coeff([y_int1 y_int2],[y_fit1 y_fit2]); % correlation coefficient
    figure(ci)
    plot(x_fit1, y_fit1, '-',x_fit2, y_fit2, '-',x,y(ci,:), 'r .', 'MarkerSize', 25)
%     title(['Gaussian Fit / R1² = ' num2str(Rsquared1) ' / R2² = ' num2str(Rsquared2)], 'FontSize', 15)
    title(['Gaussian Fit / R² = ' num2str(Rsquared_global)], 'FontSize', 15)
    ylabel('Y', 'FontSize', 14)
    xlabel('X', 'FontSize', 14)
    eqn1 = " y = "+a_sol1(ci)+ " * exp(-(x - " +b_sol1(ci)+")² / (" +c_sol1(ci)+ ")²";
    eqn2 = " y = "+a_sol2(ci)+ " * exp(-(x - " +b_sol2(ci)+")² / (" +c_sol2(ci)+ ")²";
    legend(eqn1,eqn2,'data')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function     [x_fit,y_fit,sol] = myfit(x,y)
    
    f = @(a,b,c,x) a.*exp(-(x-b).^2 / c.^2);
    obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
    sol = fminsearch(obj_fun, [max(y),max(x)/2,max(x)/6]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    x_fit = linspace(min(x),max(x),100);
    y_fit = f(a_sol, b_sol,c_sol, x_fit);
    
    
end
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

Mathieu NOE 2021-10-6

在 MATLAB Online 中打开

This is my suggestion for 2 peaks single fit

clearvars
clc
% load data
x = csvread('time.csv');
y = csvread('data_che.csv');
for ci = 1:5 % 5
    
    % curve fit using fminsearch
    [x_fit,y_fit,sol] = myfit((x),y(ci,:));
    a_sol(ci) = sol(1);
    b_sol(ci) = sol(2);
    c_sol(ci) = sol(3);
    d_sol(ci) = sol(4);
    e_sol(ci) = sol(5);
    f_sol(ci) = sol(6);
    y_int = interp1(x,y(ci,:), x_fit);
    
    Rsquared = my_Rsquared_coeff(y_int,y_fit); % correlation coefficient
    figure(ci)
    plot(x_fit, y_fit, '-',x,y(ci,:), 'r .', 'MarkerSize', 25)
    title(['Gaussian Fit / R² = ' num2str(Rsquared)], 'FontSize', 15)
    ylabel('Y', 'FontSize', 14)
    xlabel('X', 'FontSize', 14)
    eqn = " y = "+a_sol(ci)+ " * exp(-(x - " +b_sol(ci)+")² / (" +c_sol(ci)+ ")² +  "+d_sol(ci)+ " * exp(-(x - " +e_sol(ci)+")² / (" +f_sol(ci)+ ")²";
    legend(eqn,'data')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function     [x_fit,y_fit,sol] = myfit(x,y)
    
    f = @(a,b,c,d,e,f,x) a.*exp(-(x-b).^2 / c.^2) + d.*exp(-(x-e).^2 / f.^2);
    obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4), params(5), params(6),x)-y);
    
    % peaks  y values
    [PKS,LOCS]= findpeaks(y);
    sol = fminsearch(obj_fun, [PKS(1),x(LOCS(1)),x(LOCS(1))/3, PKS(2),x(LOCS(2)),x(LOCS(2))/3]); % to be updatd
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    d_sol = sol(4);
    e_sol = sol(5);
    f_sol = sol(6);
    x_fit = linspace(min(x),max(x),200);
    y_fit = f(a_sol, b_sol,c_sol,d_sol, e_sol,f_sol, x_fit);
    
end
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

Mathieu NOE 2021-10-7

在 MATLAB Online 中打开

hello

the shape of your data does not look like a periodic signal , so you can always convert your time domain signals to frequency using fft (and you can convert back to time using ifft , but the length of the fft data is equal to your time data because you need almost all the fft frequencies to get back (fit) correctly to the time data

as an example , here I have limited the fft to 3 frequencies and the fit is not very good , even if I remove the trailing zeros after t = 3500;

code :

clearvars
clc
% load data
x = csvread('time.csv');
y = csvread('data_che.csv');
ind = find(x>3500);
x(ind) = [];
y(:,ind) = [];
for ci = 1:5 % 5
    
    % curve fit using fminsearch
    [x_fit,y_fit,sol] = myfit((x),y(ci,:));
    a_sol(ci) = sol(1);
    b_sol(ci) = sol(2);
    c_sol(ci) = sol(3);
    d_sol(ci) = sol(4);
    e_sol(ci) = sol(5);
    f_sol(ci) = sol(6);
    g_sol(ci) = sol(7);
    w_sol(ci) = sol(8);
    y_int = interp1(x,y(ci,:), x_fit);
    
    Rsquared = my_Rsquared_coeff(y_int,y_fit); % correlation coefficient
    
    figure(ci)
    plot(x_fit, y_fit, '-',x,y(ci,:), 'r .', 'MarkerSize', 25)
    
    title(['Fourier Series Fit / R² = ' num2str(Rsquared)], 'FontSize', 15)
    ylabel('Y', 'FontSize', 14)
    xlabel('X', 'FontSize', 14)
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function     [x_fit,y_fit,sol] = myfit(x,y)
% y: a + b*cos(w*x) + c*cos(2*w*x) + d*cos(3*w*x) + e*sin(w*x) + f*sin(2*w*x) + g*sin(3*w*x)
f = @(a,b,c,d,e,f,g,w,x) a + b*cos(w*x) + c*cos(2*w*x) + d*cos(3*w*x) + e*sin(w*x) + f*sin(2*w*x) + g*sin(3*w*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4), params(5), params(6), params(7), params(8),x)-y);
% peaks  y values
[PKS,LOCS]= findpeaks(y);
period = x(LOCS(2)) - x(LOCS(1)); 
w_init = 2*pi*1/period;
sol = fminsearch(obj_fun, [zeros(1,7) w_init ]); % to be updatd
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
d_sol = sol(4);
e_sol = sol(5);
f_sol = sol(6);
g_sol = sol(7);
w_sol = sol(8);
x_fit = linspace(min(x),max(x),200);
y_fit = f(a_sol,b_sol,c_sol,d_sol,e_sol,f_sol,g_sol,w_sol,x_fit);
end
function Rsquared = my_Rsquared_coeff(data,data_fit)
% R2 correlation coefficient computation
% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);
% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);
% definition of the coefficient of correlation is
Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
end

Mathieu NOE 2021-10-7

在 MATLAB Online 中打开

this is a fft based solution, where you can decide to truncate the fft output to the nfft first terms

seems to work pretty well with nfft = 8 (min)

clearvars
clc
% load data
x = csvread('time.csv');
y = csvread('data_che.csv');
%fft truncation to nfft first terms
 nfft = 8;
     
for ci = 1:5 % 5
    
    % fft fit using fminsearch
    yfft = fft(y(ci,:));
    l = length(yfft);
    
    %fft truncation to nfft first terms (do not forget to generated a
    %double sided fft vector)
     db_side_fft = zeros(1,l);
     db_side_fft(1:nfft) = yfft(1:nfft);
     db_side_fft(l-nfft+2:l) = yfft(l-nfft+2:l);
     y_fit = real(ifft(db_side_fft));
    Rsquared = my_Rsquared_coeff(y(ci,:),y_fit); % correlation coefficient
    figure(ci)
    plot(x, y_fit, '-',x,y(ci,:), 'r .', 'MarkerSize', 25)
    title(['Fourier Series Fit / R² = ' num2str(Rsquared)], 'FontSize', 15)
    ylabel('Y', 'FontSize', 14)
    xlabel('X', 'FontSize', 14)
end
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
end

AS 2021-10-22

在 MATLAB Online 中打开

I want creat the matrix c from kmin to kmax.like c(2,5,1)....c(3,5,2).....c(20,5,19)....

But, using this code it is giving the final value of kmax with varying 1: 19.

So, how it will be solved. knidly, suggest me to resolve it. Thank you.

kmin=2;
kmax=20;
p=randperm(size(data,1)); % data is 3375 by 5
% generate initial centre
for i=kmin:kmax
    for i1=1:5
        for i2=1:19
        c(i,:,i2)=(data(p(i),:)); 
        end
    end
end

Mathieu NOE 2021-10-25

hello

seems to me the index i1 is not usedin the line

c(i,:,i2)=(data(p(i),:));

so It should be

c(i,i1,i2)=(data(p(i),:));

请先登录，再进行评论。

How to fit curve using for loop?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

20 个评论
显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to fit curve using for loop?

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

20 个评论 显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论