Curve fitting of two equations

Hello,

I have two equations

f=0:0.2:3;

y=teta*exp((a+1)*f);

x=exp((a-1)*f-1)/(a-1)^2;

I have a set of x vs. y

x=[36.78 ,37.53 ,38.28 ,39.06 ,39.85 ,40.65 ,41.47 ,42.3162082317748,43.17 ,44.04 ,44.93 ,45.84 ,46.76 ,47.71 ,48.67 ,49.65 ]

y=[0.01 ,0.0152 ,0.023 ,0.035,0.0536 ,0.0816 ,0.12 ,0.1891 ,0.287 ,0.438 ,0.6668 ,1.01494 ,1.544 ,2.35097424365239,3.578 ,5.445 ]

and I want to fit the data points and find out the constants teta,a.

I will appreciate any help.

Thanks

[Alex Sha]

refer to the results below:

Root of Mean Square Error (RMSE): 0.00236810622550458
Sum of Squared Residual: 0.000179453667048754
Correlation Coef. (R): 0.9999997274333
R-Square: 0.999999454866675
Adjusted R-Square: 0.999999318583346
Determination Coef. (DC): 0.999999412360353
F-Statistic: 21897434.6722157

Parameter          	Best Estimate
--------------------	-------------
teta	0.0099986478526145
a	1.100005339318

Answer 1

[John D'Errico]

I'll expand on this answer in the morning. But, yes, you can use lsqcurvefit here, or lsqnonlin.

The y equation in particular suggests logging it however. Anytime that you see large differences in scale across the range of your problem, it makes sense that taking a log will be beneficial. Essentially this presumes proportional error in the fit is appropriate.

If you don't log the y expression, then the problem is, lsqnonlin presumes additive, homoscedastic error. So it assumes that the point with y == 0.01 will have the same absolute importance as the point where y=5.445.

As well, the exponentials in there also suggest taking logs.

x = [36.78 ,37.53 ,38.28 ,39.06 ,39.85 ,40.65 ,41.47 ,42.3162082317748,43.17 ,44.04 ,44.93 ,45.84 ,46.76 ,47.71 ,48.67 ,49.65 ]

y = [0.01 ,0.0152 ,0.023 ,0.035,0.0536 ,0.0816 ,0.12 ,0.1891 ,0.287 ,0.438 ,0.6668 ,1.01494 ,1.544 ,2.35097424365239,3.578 ,5.445 ]

f=0:0.2:3;

So, you have TWO equations that will allow you to predict x and y, as a function of f, and the parameters teta and a.

y=teta*exp((a+1)*f);
x=exp((a-1)*f-1)/(a-1)^2;

Logging the y equation gives us:

log(y) = log(teta) + (a+1)*f;

If we had no more than tis equation, then we could compute teta and a directly using simple linear regression.

The x equation need not be logged, because all the values of x are essentially the same size. So additive versus proportional error is not a real factor. I'll leave it alone, although there is no problem either way.

So, now lets formulate an objective function. I'll use lsqnonlin here,as I prefer it for this specific case. It is not a strong preference though.

I'll implicitly put teta into params(1), a as params(2).

fun = @(params) [x - exp((params(2)-1)*f-1)/(params(2)-1)^2, log(y) - log(params(1)) - (params(2)+1)*f];

Now, just call lsqnonlin on fun.

You probably want to put lower bounds on teta as 0, as you don't want it trying to try out teta less than zero. In fact, I'd suggest setting the lower bound for teta as something like 1e-10, to prevent it from even trying to go negative. If it does, lsqnonlin will fail, due to the complex numbers thereby created.

Similarly, set the lower bound on a as something like -1+1e-10.

The essential idea here is to pack the x and log(y) equations together into one set, then using lsqnonlin to solve the entire set for the combined paramters teta and a.

Good starting values, just be looking at the equations MIGHT be easily gained, from the log(y) equation. Simply fit a linear polynomial to that curve.

p0 = polyfit(f,log(y),1)
p0 =
      2.102      -4.6117

That implies

a+1 = 2.102

and

log(teta) = -4.6117

solving for a and teta, we get

teta0 = exp(p0(2))
teta0 =
    0.0099353

a0 = p0(1) - 1
a0 =
        1.102

These are logical numbers. we know that a MUST be slightly smaller or larger than 1. As a test, I'll solve the problem twice, starting a from slightly above or below 1.

params = lsqnonlin(fun,[teta0,a0],[1e-10,-1+1e-10])

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the function tolerance.

<stopping criteria details>

params =
    0.0099644          1.1
>> 
>> params = lsqnonlin(fun,[teta0,0.99],[1e-10,-1+1e-10])

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.

<stopping criteria details>

params =
     0.013194      0.91282

[friet]

Wow! Thanks for you explanation. I still didnt get the first part but based on your explanation I type the following code and have some issue with the initial guess.

    clear all
  clc
  x = [36.78 ,37.53 ,38.28 ,39.06 ,39.85 ,40.65 ,41.47 ,42.31,43.17 ,44.04 ,44.93 ,45.84 ,46.76 ,47.71 ,48.67 ,49.65 ];
  y = [0.01 ,0.0152 ,0.023 ,0.035,0.0536 ,0.081 ,0.12 ,0.1891 ,0.287 ,0.438 ,0.66 ,1.01494 ,1.544 ,2.35,3.578 ,5.445 ];
  f=0:0.2:3;

fun = @(params) [x - exp((params(2)-1)*f-1)/(params(2)-1)^2, log(y) - log(params(1)) - (params(2)+1)*f];

lb = [1e-10,-1+1e-10];  % Set lower bounds 
ub = [inf,inf];         % Set upper bounds
params0 = [0.1,0.1];    % Set initial point.

options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');

params = lsqnonlin(fun,params0,lb,ub)

figure(1)
plot(x,y,'ko',exp((params(2)-1)*f-1)/(params(2)-1)^2,params(1)*exp((params(2)+1)*f),'b-')
legend('Data','Best fit')

It didnt work well :(

When i use those value it works well.

params0 = [0.1,1.1];

So my question is why the lsqnonlin is so dependent on the initial values. I was even thinking to just say 0 for both initial values and the lsqnonlin just calculate after it got the perfect a0 and teta0 that fits well.

And finally instead of using

plot(x,y,'ko',exp((params(2)-1)*f-1)/(params(2)-1)^2,params(1)*exp((params(2)+1)*f),'b-')

is there a way just to call the function and plot.

Thanks again!!

[John D'Errico]

ALL optimizers are dependent on the initial guess. ALL. Must I repeat myself? ALL OF THEM. Even "global" methods have issues.

So, think about these functions! LOOK AT THEM. After all, I did, and I saw immediately that if you are not careful at how you start them out, you MAY get trash.

So now, look at the equation for x.

x = exp((a-1)*f-1)/(a-1)^2;

What happens if a is not close to 1? After all, f varies by a fair amount, and you have an exponential function of f in there! Since x is almost constant over the entire range of f, this means that for it to be true, then (a-1) MUST be very near zero, as (a-1) is the coefficient of f.

LOOK AT YOUR EQUATIONS. UNDERSTAND THEM. They are after all, YOUR equations, not mine. You should know what is in them.

Next, you need to understand that optimizers, when started with poor starting values, will often tend to function poorly.

I've said this many times before. Think about an optimizer as a blind person, set down on the surface of the earth, and tasked with finding the lowest spot. All they can do is to tap their cane and move in a direction they think will help, based on the local environment, and what they have seen before. (I'll be nice, and allow the person to carry scuba gear.) If you start this fellow in a bad place, will you seriously expect success?

Finally, as for plotting the curves, you have a very simple solution. What is wrong with it?

[friet]

Thanks for your comment. Yes! You are right, The equations are understandable, unfortunately, I am dealing with very complected that those I brought for discussion.

Actually after I read more about this "lsqnonlin" solver in matlab I found out if you set the solver to follow the so called levenberg-marquardt algorithm you can enter a wide range of initial solution. In other words the solver will converge regardless of what initial value you select.

options.Algorithm = 'levenberg-marquardt';

params = lsqnonlin(fun,params0,lb,ub,options)

So i am all set with this

Thanks

Finally i couldn't figure out the very simple solution for plotting from function.

Many Thanks

Answer 2

[Walter Roberson]

You do not have two equations: each value of f generates another pair of equations.

It appears there might be two solutions for teta

teta1 = mean( exp((1/2)*(-f^2+2*x*ln(y)+sqrt(f^4-4*f^2*x))/x) )
teta2 = mean( exp((1/2)*(-f^2+2*x*ln(y)-sqrt(f^4-4*f^2*x))/x) )

a = mean( (ln(y)-ln(teta)-f)/f )   for teta = teta1 and teta = teta2

No proof offered at this time; my eyes are giving me trouble.

[friet]

Thanks for your answer. But I do have two equations. I explain the question more. I think it has to do sth with "lsqcurvefit"