Linearization using the power function
15 views (last 30 days)
Show older comments
Hello !
I got these data
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]
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?
0 Comments
Answers (1)
Walter Roberson
on 25 May 2023
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)
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)
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
0 Comments
See Also
Categories
Find more on Linear and Nonlinear Regression in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!