nonlinear curve fitting a function on data

1 visualización (últimos 30 días)
Moj
Moj el 9 de Jun. de 2019
Editada: John D'Errico el 28 de Feb. de 2020
Hi,
I have some experiment data. Hereby, I need to fit the following function to determine one of the variable. A Levenberg–Marquardt least-squares algorithm was used in this procedure.
My experiment data:
beta = 1.135e-4;
sin(theta) = [-0.81704 -0.67649 -0.83137 -0.73468 -0.66744 -0.43602 0.45368 0.75802 0.96705 0.99717 ];
x = [72.01 59.99 51.13 45.53 36.15 31.66 30.16 29.01 25.62 23.47 ];
And function which is needing to be fitted:
sin(theta) = -1+2*sqrt(alpha/x)*exp(-beta*(x-alpha)^2)
Is there any suggestion to determine alpha variable?
Thanks in advance

Respuesta aceptada

John D'Errico
John D'Errico el 9 de Jun. de 2019
Editada: John D'Errico el 9 de Jun. de 2019
sin(theta) is not a variable name.
beta = 1.135e-4;
sintheta = [-0.81704 -0.67649 -0.83137 -0.73468 -0.66744 -0.43602 0.45368 0.75802 0.96705 0.99717 ];
x = [72.01 59.99 51.13 45.53 36.15 31.66 30.16 29.01 25.62 23.47 ];
Before you do ANYTHING, PLOT your data. Ask, if it is consistent with your model. The answer here, is not really.
plot(x,sintheta,'o')
grid on
axis([20 80 -1 1])
yline(-0.8);
I drew a line at y = -0.8, which is where I would suggest the baseline for your curve should roughly lie at.
That is, what is the asymptotic value of your function as x approaches infinity? You have written it as:
sin(theta) = -1+2*sqrt(alpha/x)*exp(-beta*(x-alpha)^2)
As x gets large, this approaches -1. Yet, clearly, your data does not get anywhere near -1, and it never will. Your data approaches a value closer to -0.8, MAYBE -0.85. So your data is inconsistent with your model. That you want to call the depedent variable sin(theta) does not matter. It simply will never get anywhere near -1.
Lets look at how you might fit this with the curve fitting toolbox. remember, that you will NEED good starting values for alpha and beta.
mdl = fittype('-1+2*sqrt(alpha./x).*exp(-beta*(x-alpha).^2)','indep','x');
fittedmdl = fit(x',sintheta',mdl,'startpoint',[20,0.01])
fittedmdl =
General model:
fittedmdl(x) = -1+2*sqrt(alpha./x).*exp(-beta*(x-alpha).^2)
Coefficients (with 95% confidence bounds):
alpha = 25.82 (23.66, 27.99)
beta = 0.01954 (0.0005392, 0.03854)
plot(fittedmdl)
hold on
plot(x,sintheta,'o')
yline(-.8);
Note that I estimated beta there, too, as the value for beta that you gave is also inconsistent with the data. If I try to force the model to use beta =1.135e-4;, I get garbage for a fit, far worse than what I got for the model I chose.
Now, look carefully at the curve that came out. It does not look at all like your data. The baseline is wrong. The shape of the curve in the peak is wrong.
Why do I say those shapes are "wrong"? The peak of the curve is fairly wide. It rolls over nice and slowly. And it never even manages to start coming back down, if I look at your data. But the transition region where it drops down from 1 towards -1, is VERY sharp. That means you need a relatively large value for beta, NOT one as small as you seem to think it should be. Otherwise, you cannot get a sharp transition. But with a sharp transition there, that also means the peak wants to drop down for small x, below the value of alpha. So your data is simply not consistent with that model.
I honestly don't know where you got that model. But it simply will never look like that data. Sorry. But trying to fit a square peg into a round hole tends only to damage either the peg or the hole.
(In fact, if I try to estimate the baseline from your data, it comes out to approximately -0.79, so my visual guess was pretty good.)
  2 comentarios
Moj
Moj el 28 de Feb. de 2020
Could you also please put any comments regarding how is possible to calculate fitting error?`
John D'Errico
John D'Errico el 28 de Feb. de 2020
Editada: John D'Errico el 28 de Feb. de 2020
What is the problem?
First, what is "error"? Logically, it would be the difference between the measured data values, and the predicted values. Then you might compute some measure of the errors for the entire data set. That might mean to compute the standard deviation of the errors.
If you are using the curve fitting toolbox, this is all easy.
In the example I gave in my answer, you can evaluate any point or list of points with the simple exppresion:
ypred = fittedmdl(x);
yerr = y - ypred;
errstd = std(yerr);
So the list of errors are now the vector yerr. You can now trivially compute the standard deivation of the errors. Perhaps slightly better is to compute an RMSE, where you divide the sum of the squares of the errors by the number of data points minus the number pf parameters estimated. But WTP? In the model above, there were 2 parameters estimated, thus alpha and beta.
RMSE = sqrt(sum(yerr.^2)/(numel(y) - 2));

Iniciar sesión para comentar.

Más respuestas (1)

Matt J
Matt J el 9 de Jun. de 2019
Editada: Matt J el 10 de Jun. de 2019
Are you sure the model isn't supposed to be a Gaussian+constant? It gives a much better fit. I used this FEX file to do the fit, but you could use the Curve Fitting toolbox as well, if you have it.
LB = [];
UB = {[],[],xdata(end)} ; %center the Gaussian left of the first data point
params0={-1,2,30,[]};
[params] = gaussfitn(xdata,ydata,params0,LB,UB);
[D,A,mu,sig]=deal(params{:})
fun=@(x) D + A*exp( -0.5 * (x-mu).' * inv(sig) *(x-mu) );
xc=linspace(xdata(1),xdata(end),1000);
plot(xdata,ydata,'o',xc,arrayfun(fun,xc),'-')

Categorías

Más información sobre DSP Algorithm Acceleration en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by