Fitting gaussian exponential to logscale

Hello!

I have been trying to fit a log scale plot to an exponential function given below:

fun=@(t,x)2*(t(1))^2.*(1-exp(-((x./t(2)).^(2*t(3)))));

But I am unable to get a fit that is even close to the data plot. this is what I tried to do:

function fun = myfun(t,x,y)
t(1) = t(1)*1e-9;
t(2) = t(2)*1e-9;
t(3) = t(3)*1e-9;
xx=log(x)
yy=log(y)
fun=@(t,x)2*(t(1))^2.*(1-exp(-((x./t(2)).^(2*t(3)))));;

t=lsqcurvefit(fun,t,x,y)

plot(xx,yy,'ko',xx,fun(t,xx),'b-')
end

and then caling the function as:

t=[7.25553e-10 2.2790e-9 2.27908e-9]

myfun(t,x,y)

But this is what I get instead:

The x and y data is attached as .txt

I tried playing with different inital values but none of them gives me close to a good fit. I dont understand what is wrong here.

Could someone please help me out ?

I also do not have access to the curve fitting tool or any other tool in matlab just so you know. Any help would be really appreciated!

Thank you!

Answer 1

[Alan Stevens]

The logs of the y values are negative, whereas your function will return positive values only! You are probably better off trying the following

% Scale data
xx = 10^7*x;
yy = 10^18*y;

fun=@(t,x) 2*(t(1))^2.*(1-exp(-((x./t(2)).^(2*t(3)))));

t0 = [1 1 0.5];
t = fminsearch(@(t) fcn(t,xx,yy,fun),t0);
disp(t)
Y = fun(t,xx);
% Rescale curve fit
yfit = Y*10^-18;

plot(x,y,'o',x,yfit),grid
xlabel('x [m]'),ylabel('y [m^2]')
legend('data','curve fit')

function F = fcn(t,xx,yy,fun)

         Y = fun(t,xx);
         F = norm(Y-yy);

end

which results in the following fit

[sandy]

Thank you! this works well!

Just a question, do you know how can I find the fit parameter error range in +/- ?

[Alan Stevens]

I don't think you can get that automatically from fminsearch. However, it's possible that taking the same approach with lsqcurvefit, or lsqnonlin, the error range might be returned automatically (I can't tell as I don't have the Optimization toolbox).

Answer 2

[Mathieu NOE]

hello

this a code that seems to work

hope it helps

I prefered to simplify as much the work of the optimizer , so first convert x and y in log scale and shift origin so curve start at x = y = 0

clc
clear all
close all

 T = readtable('trial.txt');
 C = table2array(T);
 x = C(2:end,1);
 y = C(2:end,2);
 
 % convert x and y in log scale
 lx = log10(x);
 ly = log10(y);

 % center origine so lx and ly first point = 0
 mlx = min(lx);
 lxc = lx - mlx;
 mly = min(ly);
 lyc = ly - mly;
 
% resample for higher points density in the first half of lx
% seems to help get better results
 lxxc = linspace(min(lxc),max(lxc),length(lxc));
 lyyc = interp1(lxc,lyc,lxxc);
 

 plot(lxc, lyc,'+b',lxxc, lyyc,'+r');

f = @(a,b,x) a.*(1-exp(b.*x));
obj_fun = @(params) norm(f(params(1),params(2),lxxc)-lyyc);
sol = fminsearch(obj_fun, [1,1]);
a_sol = sol(1);
b_sol = sol(2);

figure;
plot(lxc, lyc, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
plot(lxxc, f(a_sol, b_sol, lxxc), '-');grid

% finally 
% lyc = a.*(1-exp(b.*lxc)); % converted back to : 
% ly = mly + a.*(1-exp(b.*(lx - mlx)));
% and in linear scale : 
y_fit = 10.^(mly + a_sol.*(1-exp(b_sol.*(log10(x) - mlx))));
figure;
plot(x, y, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
plot(x, y_fit, '-');grid

[sandy]

Thank you for the code. I think it works well for a fit but I believe you use a different function than the one I mentioned in the post ?

Also, I am not quite clear on the following part, is it to improve the inital parameter guess for y_fit function ?

f = @(a,b,x) a.*(1-exp(b.*x));
obj_fun = @(params) norm(f(params(1),params(2),lxxc)-lyyc);
sol = fminsearch(obj_fun, [1,1]);
a_sol = sol(1);
b_sol = sol(2);

[Mathieu NOE]

hello

sure , there are more than one solution to the problem, but as I am lazzy, I stick to the solutions I know that worked in the past (un less I'm forced to get out of my comfort zone !)

Also, I am not quite clear on the following part, is it to improve the inital parameter guess for y_fit function ?

Which part is not clear ?

tx

[sandy]

The part in the snippet of code that I mentioned in the previous comment, I wanted to be clear, is the use of a function

f = @(a,b,x) a.*(1-exp(b.*x));

just to improve the parameter for yfit function?

[Mathieu NOE]

no, it's basically where you define the model fit equation

the fminsearch will use that funtion in it's evaluation

for more details , look at : help fminsearch

tx