# Fit power function to data (optimization)

17 views (last 30 days)
Alfredo Scigliani on 7 Dec 2021
Commented: Star Strider on 7 Dec 2021
I am trying to fit a general power function ydata=a*(xdata)^b (where a and b are fitting parameters, xdata and ydata are the given data).
I am trying to use lsqcurvefit to optimize for the parameters a and b but I don't believe I am getting the best fit possible. I believe that the fitting cout be way better, maybe lsqcurvefit is not the best idea for fitting a power function.
Matlab Script:
--------------------------------------------------------------------------------------------------------------------------------------------------------------
clear; clc; close all;
xdata = [0.02 0.05 0.1 0.2 0.3 0.4 0.5];
ydata = [1.7935e-006 5.6359e-006 10.8276e-006 38.1193e-006 60.4152e-006 65.8997e-006 62.0615e-006];
fun = @(x,xdata) x(1)*xdata.^x(2);
x0 = [1,0.5]; %initial guess
x = lsqcurvefit(fun,x0,xdata,ydata)
loglog(xdata,ydata,'ro', xdata, fun(x,xdata),'r-')
-------------------------------------------------------------------------------------------------------------------------------------------------------------- As you can see, the fitted line could be way better and pass through more data .
If you have any suggestion in how I can improve my code to obtain a better fit, I am open to it and I appreciate it.

John D'Errico on 7 Dec 2021
Your problem is, you are trying to use LEAST squares. But look carefully at your data.
xdata = [0.02 0.05 0.1 0.2 0.3 0.4 0.5];
ydata = [1.7935e-006 5.6359e-006 10.8276e-006 38.1193e-006 60.4152e-006 65.8997e-006 62.0615e-006];
Do you see that some of your data is nearly 2 powers of 10 larger then the rest? You plot it on a log scale, so you do't really think of it. It looks like just a straight line on those axes.
So what happens is error in some of your data points are WAY more important than are errors in the others. This is a case where you do NOT want to use a simple least squares to estimate those coefficients. You might decide to use a weighted least squares. A simple alternative is to use a straight line fit, to the log of your data.
Why is that a good idea? Because in the case you have, you really do not want to use an error structure that is additive, with a uniform variance on all data points. Instead, you want to use a proportional error structure, s multiplicative errors. What happens when you take the log? It turns what you have, thus a lognormally distributed error, into normally distributed, additive error.
xlog = log(xdata);
ylog = log(ydata);
P1 = fit(xlog',ylog','poly1')
P1 =
Linear model Poly1: P1(x) = p1*x + p2 Coefficients (with 95% confidence bounds): p1 = 1.189 (0.9811, 1.397) p2 = -8.535 (-9, -8.071)
plot(P1,xlog,ylog,'o') I used fit there, because it makes the plotting really easy.
But you should see how well it worked here, giving you the model you were hoping to find.
Now we need to transform the model back into a and b.
b = P1.p1
b = 1.1890
a = exp(P1.p2)
a = 1.9645e-04
The distinction I pointed out here is an important one.
Star Strider on 7 Dec 2021
Alex Sha uses a completely different (and expensive) optimisation package, not MATLAB.
Ths point is that a power law fit does not work with the available data. It’s simply not an appropriate model.

### More Answers (1)

Star Strider on 7 Dec 2021
The problem is the model. It may not be appropriate for these data.
A logistic curve fits much better —
xdata = [0.02 0.05 0.1 0.2 0.3 0.4 0.5];
ydata = [1.7935e-006 5.6359e-006 10.8276e-006 38.1193e-006 60.4152e-006 65.8997e-006 62.0615e-006];
% fun = @(x,xdata) x(1)*xdata.^x(2);
fun = @(b,x) b(1)./(1+exp(b(2).*(x-b(3))));
% x0 = [1,0.5]; %initial guess
x0 = rand(1,3).*[1E-5 -15 median(xdata)];
x = lsqcurvefit(fun,x0,xdata,ydata)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
x = 1×3
0.0001 -23.1434 0.1804
figure
plot(xdata,ydata,'ro', xdata, fun(x,xdata),'r-') figure
loglog(xdata,ydata,'ro', xdata, fun(x,xdata),'r-') .

R2019b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!