Linearization using the power function

Hello !

I got these data

x = 0:0.2:1.8

x = 1×10

0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000

y = [1.3 2.1 3.0 5.2 8.4 13.5 22 33.5 53 85.4] 

y = 1×10

1.3000 2.1000 3.0000 5.2000 8.4000 13.5000 22.0000 33.5000 53.0000 85.4000

and I'm being asked to find the best fit function for data. So I decided to use power function for the linerization but when I'm calculating the a1 and a0 ( y = a1*x ^(a0) ) i dont have any result (a1 =a0 = NaN) .

clear;
close all;
clc;
format short;

x = 0:0.2:1.8;
y = [1.3 2.1 3.0 5.2 8.4 13.5 22 33.5 53 85.4];

figure(1)
subplot(1,2,1)
plot(x, y, 'o')
title('Original Data')

% Linearization using the power function
xx = log10(x);
yy = log10(y);

subplot(1,2,2)
plot(xx, yy, 'o')
title('Power function') 

[a1, a0] = LinearRegression(xx,yy);

%----------------------------------------------------------------------------------------------------------------------------------------
function [a1, a0] = LinearRegression(x,y)
n = length(x); 
Sumx = sum(x); Sumy = sum(y); Sumxx = sum(x.*x); Sumxy = sum(x.*y);
den = n*Sumxx - Sumx^2;
a1 = (n*Sumxy - Sumx*Sumy)/den; a0 = (Sumxx*Sumy - Sumxy*Sumx)/den;
% Plot the data and the line fit
l = zeros(n,1);  % Pre-allocate
for i = 1:n
    l(i) = a1*x(i) + a0;   % Calculate n points on the line
end
plot(x,y,'o')
hold on
plot(x,l)
end

Can you help me with this please?

Answer 1

[Walter Roberson]

x = 0:0.2:1.8;
y = [1.3 2.1 3.0 5.2 8.4 13.5 22 33.5 53 85.4];

x starts with zero

xx = log10(x);

log10(0) is

You cannot do meaningful linear regression with terms that include infinity.

You are aiming for a formula of y = c * 10^x . If you calculate for x == 0 then 10^0 == 1 so y(0) = c and you can then calculate

c = y(1);
yprime = y ./ c;
p = polyfit(log(x(2:end)), yprime(2:end), 1)

p = 1×2

22.6779 23.5658

but if y = c * 10^x were true then p(2) would have to equal log(1) == 0 after we divided y by y(0).

We can conclude from this that you cannot fit the data to the kind of curve you want.

How about:

xL = log10(x(2:end)); yL = log10(y(2:end));
p1 = polyfit(xL, yL, 1);
p2 = polyfit(xL, yL, 2)

p2 = 1×3

1.8501 2.4612 1.1487

yproj1 = 10.^polyval(p1, log10(x));
yproj2 = 10.^polyval(p2, log10(x));
plot(x, y, 'ko', 'displayname', 'original points');
hold on
plot(x, yproj1, 'rv', 'displayname', 'log10 linear');
plot(x, yproj2, 'b-.', 'displayname', 'log10 quadratic');

p4 = polyfit(x, y, 4);
y4 = polyval(p4, x, y);
plot(x, y, '--g+', 'displayname', 'quartic polynomial')
legend show

The red triangles, log10 linear, is the fit you were trying to use. Using a quadratic in log10, or using a plain quartic in linear space, gives much better results