7 views (last 30 days)

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

John D'Errico
on 9 Jun 2019

Edited: John D'Errico
on 9 Jun 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.)

John D'Errico
on 28 Feb 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));

Matt J
on 9 Jun 2019

Edited: Matt J
on 10 Jun 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),'-')

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.