nonlinear curve fitting a function on data

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

Answer 1

[John D'Errico]

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.)

[Moj]

Could you also please put any comments regarding how is possible to calculate fitting error?`

[John D'Errico]

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

Answer 2

[Matt J]

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),'-')