error calculation perfect fit

Hi guys, so I have experimental data in column matrix y. And I want to use the equation below to find x as a function of y. This means have constant b adjusted from 0 to 2 and hence want to obtain the ideal b value.

I'm thinking of using a double for loop however would anyone know how to calculate for error and set a desirable tolerance. Or tips on how to solve for it effectively?

Thanks in advance. The equation is:

y= (1-0.05*x^(b))/(1-0.03*x^b)^2

[Star Strider]

See also for context: fitting data with equation.

Answer 1

[Star Strider]

I would use nlinfit and nlparci to estimate ‘b’, and use a few values of ‘x’ and ‘y’ to provide an initial estimate, with:

b = log(20*(1-1./y))./log(x);

Noting the assumptions and restrictions that I mentioned earlier.

[sarah]

I'll try this method and will update you if it works or not. Thanks

[Star Strider]

My pleasure.

Please be mindful of the restrictions and assumptions.

[sarah]

Okay so now I have experimental data x and y which follows an exponential trend. i.e. as x increases y decreases.

And i have two beta variables to predict from my model b(1) and b(2). It doesn't have noise hence do you think it will give an accurate answer for b(1) and b(2)?

Just found a way to link my experimental data x and y (found new equation).

[Star Strider]

Without noise, it should give reasonably accurate estimates of both parameters. However, the model you posted has only one ‘b’ parameter. If you want to get accurate estimates for both (for instance in chemical or biochemical kinetics or some such), you have to fit and estimate both in the same nonlinear fit. Fitting each separately is known as ‘curve stripping’. It has no theoretical foundation, and leads to inaccurate parameter estimates, even with very low-noise data.

I have no idea what your data are or what you’re doing, so I can only guess. I would like to be sure you are estimating your parameters correctly, and so will get usable — and publishable — parameters from your data, if that’s what you intend to do.

[sarah]

Yeah I incorporated another equation and now have two b(1) and b(2) variables. I'll use it and see if I get reasonable estimates and let you know if I have difficulty. Thanks a lot.

[Star Strider]

I would use nlinfit to estimate the parameters with your original equations, and not the transformed ones.

[sarah]

with the nlinfit, i still want to predict b(1) and b(2) in terms of x and y however my model has to be in terms of y i.e. it be y=2*x*y? Reason being my equation is very complex and I can't rearrange for y but i still want to use nlinfit.

Any suggestions?

[Star Strider]

I would have to see your model and some of your data to give you an appropriate reply.

You might have to abandon the curve-fitting approach and go with an optimisation routine like fsolve instead, if you cannot separate your variables to different sides of your model equation.

In your example, you would use fsolve with your objective function being:

f = @(b) b(1).*x.*y - y;

in this instance estimating your one parameter.

[sarah]

                                                                                  modelfun = @(b,X)(1-0.05*(((((b(1)*(0.00000382).^(2-b(2)))./(29.6*X)).^(1./(4-b(2))))/(0.00000382)).^(3-b(2))))./(1-0.05*((((((b(1)*(0.00000382).^(2-b(2)))./(39.6*X)).^(1/(4-b(2))))./(0.00000382)).^(3-b(2)))./0.5714)).^2;

X =

    9.9999
   10.4711
   10.9648

>> Y

Y =

    6.2016
    5.9767
    5.4755

   beta0=[0.000005048;2];

[sarah]

yes model is quite complex however i'm using few data points for now and for now i'm assuming y can be rearranged however estimates for b(1) and b(2) are no way near my initial guesses; I'll check later. Can fsolve solve for more then one variable? Sorry about the messy equation its quite complex.

[Star Strider]

The fsolve function can solve for several parameters (here ‘b’) with functions of the form:

F(b,x,y) = 0

so if you can get your equations into that form (so that it or they are all equal to zero), then you can use it to estimate your parameters. That would be appropriate when you cannot separate your expression in terms of one dependent variable as a function of one or more independent variables.

I still do not understand your essential model or what you are doing with it.

[sarah]

Thanks.

Ok so I'll explain it to you. I have a complex model equation which I want to use to predict two constants b(1) and b(2). I know there typical ranges. I have experimental data X and Y column vectors which I would like my equation model to fit on hence get an estimate for the coefficients b(1) and b(2).

nlinfit gave me values which were completely different to what I provided. I'm assuming its a typo in the equation.

(1) complex model can either have Y on one side or both sides (based on scientific assumptions) So I stared with Y on one side but it still doesn't give me accurate results.

I think I'll try nlinfit first and see if I'm making typos. Then if doesn't work I'll try fsolve. However I have two coefficients which I want to estimate and don't think fsolve can estimate two coefficients together for two column vectors X and Y (experimental data). I simply want my model equation to fit the data I have provided and hence get two coefficients. Sounds easy but it's difficult. Thanks for your patience

[Star Strider]

The fsolve function can solve one or a system of equations for one or more parameters. It can definitely do what you want, if you cannot separate your equation into a dependent variable as a function of your independent variables.

Remember to vectorise your function to do element-wise array operations unless you want to do matrix operations. You seem to have done that, but it is always the first thing to check for. (When in doubt, vectorise all operations if you intend element-wise calculations.)

[sarah]

Thanks. Hence you mean convert all my column matrices to vectors im assuming to do element to element calculations. I'll give that a go, fingers crossed it works.

[Star Strider]

No, I did not mean that you convert all your column matrices to vectors. Keep them as matrices, and address the columns as necessary. That’s easy enough in MATLAB.

By ‘vectorise’, I intended that you use array operations. (You don’t seem to be doing matrix calculations, but element-wise calculations.) See the documentation for Array vs. Matrix Operations for a full discussion of the differences.

[sarah]

Thank you. I'm getting this error for some strange reason. Looking into choosing another optimum solver now

                                                  First-Order                    Norm of 
 Iteration  Func-count    Residual       optimality      Lambda           step
     0           3         73.3319        4.54e+07         0.01
     1           6         17.9357         5.6e+05        0.001         1.0877
     2           9         2.28484            23.7       0.0001          4.576
     3          21         2.27261            32.7       100000    0.000236578
     4          25         2.27045            33.7        1e+06    3.26701e-05
     5          28         2.26815            34.6       100000     3.3644e-05
     6          32         2.26574            35.4        1e+06    3.45362e-05
     7          36         2.26549            35.5        1e+07    3.53694e-06

             No solution found.

               fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance, but the vector of function values
is not near zero as measured by the default value of the function tolerance.

criteria details

b =

   4.1255 - 3.2581i
  -0.0003 - 0.0001i

[Star Strider]

My pleasure.

I am not certain what the problem is. Your equation has to be in the form such that it equals zero. Since you mentioned that you cannot solve for one variable in terms of the other, and so cannot use nlinfit, but had terms of both variables and parameters on both sides of the equation, I suggested that you rearrange your equation so that it equals zero, and then use fsolve.

I know that this is not your function, but the idea would be that if you had for instance:

sin(x+y+b(1)) = cos(x*y*b(2))

you would rearrange it so that:

sin(x+y+b(1)) - cos(x*y*b(2)) = 0

then give it to fsolve in the form:

f = @(b,x,y) sin(x+y+b(1)) - cos(x.*y.*b(2));

[B,fval] = fsolve(@(b) f(b,x,y), b0);

(Include an options structure in the fsolve call if you want one.)

I am not certain what expressions you began with, so I cannot help you put your equation (or system of equations) in a form that fsolve can estimate, but this is the essence of the idea. (You mentioned that there could be errors in the equation you posted earlier, so I have not looked at it closely.)

[sarah]

Yes i understood and hence rearranged my function in the form of F(x,y,b)=0. Are you sure it accepts data input as column matrix?

hmm i didn't define my function as myfun = @(b,X,Y) but did it as @(b) even though i included the column vectors X,Y in the equation.

So solving for it now would be [b,fval,exitflag] = fsolve(myfun,b0,options) but now its saying not enough input arguments in original equation lol

Thankyou anyways. I'll try looking at the errors and no my equation is correct, but i'll just have a look to see what's wrong. I doubt it is the algorithm, it does say fsolve cannot handle a non-square systems which is strange.

[Star Strider]

If you have one vector each for ‘x’ and ‘y’, then this syntax should work:

[B,fval] = fsolve(@(b) myfun(b,x,y), b0);

[sarah]

its working now. Thanks so much Star strider (don't know your real name). It's giving more realistic values closer to ones which i give however it is not converging to a solution i.e. fsolve stopped because the relative size of the current step is less than the default value of the step size tolerance. Thanks anyways :)

[Star Strider]

My pleasure. I’m glad it’s working.

You can set the step size tolerance and other performance characteristics to different values by creating an options structure with the optimoptions function. Pass the structure to your fsolve call as the third parameter.

For instance, your code might change to:

opts = optimoptions('fsolve', 'TolFun', 1E-8, 'TolX', 1E-8);
[B,fval] = fsolve(@(b) myfun(b,x,y), b0, opts);

and whatever other options you want to change, although the default value for 'FinDiffRelStep' is likely as low as it should go. Experiment with the others to get the result you want.

[sarah]

thank you so much. Its coverging now :) and I didn't even have to change the tolerances lol. Your awesome keep it up! Wish I could be as expert as you :)

[Star Strider]

My pleasure, as always!

One very good way to learn MATLAB is to hang out here on Answers when you have the time. That’s how I learned a lot of what I might otherwise not have learned had I simply solved my own research problems. Read the Questions and Answers on problems that interest you, and if you believe you can provide an Answer (or a better Answer than any posted) that solves the problem, write and test the code for your Answer, and post it. Over some months, you could certainly earn the 3000 Reputation points necessary to become Editor, and more.

[sarah]

Sure.

Also i wanted to ask you regarding residuals. Residuals for some selective data is very small e-8 and the system solves. However for say the whole range i.e. 70 data points the fit isn't good i.e. residuals are large 1000. However my b(1) and b(2) are satisfactory. It doesn't mean my model is wrong right i.e. my coding. I'm assuming that experimental data at certain range is wrong.

[Star Strider]

It’s difficult to determine the reasons the residuals are very low in one region and high in another. My guess is that either (1) your model does not completely describe the process that created your data, or (2) some regions simply have more noise than others. This may not be a problem if your parameter estimates are acceptable and your model provides a good fit in your area of interest.

Answer 2

[John D'Errico]

How many points do you have that you expect a perfect fit? (Hint: with only one parameter to solve for, only ONE point will give you a perfect fit in general. With more points than that, you won't get a perfect fit.)

Anyway, with the model you have chosen, is there some good reason why you choose to define it as

y= (1-0.05*x^(b))/(1-0.05*x^b)^2

instead of the algebraically identical

y= 1./(1-0.05*x.^b)

(Note my use of ./ and .^. This is important.)

And with that model, I seriously doubt you will get a very good fit, and I have not even seen your data. That model has very little capability of adjustment.

[sarah]

my y data ranges from 1-6 but yeah i have alot of points.

so say if my experimental data is y. And say if I want to solve for b ie. the most ideal value of b. So there is no way i could do that? I think even a point of intersection which shows one value of b is good for me. As long as i get a good estimate of b - i was initially thinking of a surface plot. I do know that x ranges from 2-40. haha