function to produce line graph

Hello,

I am trying to produce this kind of line graph (sample.png). What kind of function should i use for this? I tried checking on tools > basic fitting and the 7th degree polynomial line produce almost similar line graph (plot.png), but with a note 'badly conditioned'. I also tried polyfit function, but it gave me a big curve, very different from the sample.

Many thanks for your help!

[Walter Roberson]

Be sure to use polyfit with the extra two outputs, to tell it to center and scale the data before doing the fitting. And be sure to pass those two extra values to polyval() so it can interpolate properly.

[ahamahoru]

Thank you for the suggestion, I tried to add the two outputs but those gave me a slight difference between that and plot.png result. I found 'fit' function, however, and by using extra 'power1' input, it gave me the exact same line graph as the sample along with the coefficients i needed.

Again, thanks for giving me the idea!

Answer 1

[David Wilson]

I would strongly suggest you don't try a polynomial fit here. From the context (which you neglected to tell us), I suspect that you havge a good idea what an underlying function should be. Looking at the real plot (sample.PNG), I'm expecting a 1/x or perhaps1/x^2 sort of function.

There are a number of ways to fit such functions, but a flexible way is to use the nonlinear curvefitting routines from the optimisation toolbox. That way, you can adjust the parameterisation as you see fit.

Here's my attempt at your data:

(Since you neglected to give us the data, I've taken rough estimates from your sample plot.)

x = [500, 1000, 1500, 1500, 2500, 3000, 5000, 6000, 11000, 14000, 15000, 16000]; 
y = [275, 100,  140,  40,   10,   16,   2,    4,    5,     1.5,   0.2, 0.2]; 

plot(x,y,'r>')

It's hard to see, but I 'm guessing the actual "y" data is never negative.

For my first go, I'm going to try and fit the function

This gives me two parameters, and I'm hoping p_2 will be very small (and positive). Note: there are many other functions I could have used, and you will see acouple below.)

%% Now use lscurvefit 
f = @(p,x) p(1)./(x) + p(2);  
p0 = [1e4, 0]
p = lsqcurvefit(f, p0, x,y); 

xi = linspace(1, 20000, 1e3)';
yi = f(p,xi); 

plot(x,y,'r>',xi,yi)
yline(0,'b--')
ylim([-10 300])

Actually this isn't too good. The fitted function goes negative. (If that's OK, then stop readng now.)

We could enforce p_2 to be nonnegative with some lower bounds:

f = @(p,x) p(1)./(x) + p(2);  
p0 = [1e5, 0]; 
UB = [inf 10]; % upper bounds
LB = [0   0];  % lower bounds 
p = lsqcurvefit(f, p0, x,y, ...
                LB, UB); 

xi = linspace(1, 20000, 1e3)';
yi = f(p,xi); 

plot(x,y,'r>',xi,yi); yline(0,'b--'); ylim([-10 300])

But I don't really like the resulting fit. Why not try something else, say

f = @(p,x) p(1)./(x.^2) + p(2);  
p0 = [1e5, 0]; 
UB = [1e10 1]; 
LB = [0   0]; 
p = lsqcurvefit(f, p0, x,y, ...
                LB, UB); 

xi = linspace(1, 20000, 1e3)';
yi = f(p,xi); 
plot(x,y,'r>',xi,yi); yline(0,'b--'); ylim([-10 300])

which gives a reasonable fit.

Whether that's acceptable, only you, or a domain specialist for your application can really say. I'd say it's OK, and certainly illustrates better behaviour than a polynomial.

Note that I have not taken any particular care in scaling or normalising; (so I live dangerously in these unpresented times). To be prudent, one should scale and range this problem.

[ahamahoru]

Thank you so much for your efforts to answer this question! I really appreciate that! However I found the 'fit' function with 'power1' input to be effective, and i got the coefficients that i needed as well as the R square value. But your explanation is interesting and i might need to use that in the future. For that I have some questions, could you please show me the explanation link (or just a keyword to search for) for this function f = @(p,x) p(1)./(x.^2) + p(2)? was the purpose of @(p,x) to create a new fitting equation? Never knew such thing exists until now. Also, how can i calculate R square from that?

Thank you!

[David Wilson]

Take a look at anonymous functions,

https://au.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html

To calculate R^2, check out

https://en.wikipedia.org/wiki/Coefficient_of_determination

yhat = f(p,x); % predicted values from the model  
SSE = sum((y-yhat).^2)
SST = sum((y-mean(y)).^2)
R2 = 1 - SSE/SST