How to fit curve using for loop?

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

end

[Mathieu NOE]

hi

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

tx

Answer 1

[Mathieu NOE]

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

[AS]

Thank you for the solution. I am attaching 5 curves with time. These are fitting with gauss4 with good R-square: 0.9979 value in curvefitting apps. But, in this code its not not fitted properly. I want to fit these curve in loop.

Kindly suggest me to resolve this. Thank you.

[Mathieu NOE]

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

[AS]

Thank you for the solution.For each curves hoe to store the coefficients values.

[Mathieu NOE]

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

[AS]

Thank you for the solution.

I have one more query.

[Mathieu NOE]

yes ?

[AS]

Without separating the curves into first and secon peak, is it to possible to fit the curve.

[AS]

can it be possible to fit all curves without separating them by using Fourier series?

[Mathieu NOE]

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

[AS]

Thank you for the solutions.

[Mathieu NOE]

My pleasure

[AS]

[AS]

Ig I got any doubt then i will agian ask you.

[Mathieu NOE]

ok

[AS]

Is it possible to these curves using fourier series.

[Mathieu NOE]

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]

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]

Thank you.

[AS]

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]

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