How do create a loglog line of best fit?

So far I've plotted my data and found that a loglog plot gives the most linear result. The line of best fit, however, isn't linear. How do i fix it? It also displays "negative numbers ignored" which might tie into it.

%%HK
close all % closes all windows
clear all % clae all variables
clc % clear command windows
beep off % turn beep noise off
format compact % tight spacing 

%%Data
temp = [100,200,400,600,800,1000,1200,1400,1600]
density = [3.5,1.7,0.85,0.6,0.45,0.35,0.3,0.25,0.2]

%%Plot
loglog(temp, density, '*')
xlabel('Temperature T,K')
ylabel('Density D, kg/m^3')
title('Temperature vs. Density')
grid on

%%Equation constants
Const = polyfit(temp,density, 1)
m = Const(1)
k = Const(2)

%%Bestfit Line
hold on
YBL = m*temp+k
loglog(temp,YBL)

Answer 1

[Star Strider]

Try this:

Const = polyfit(log(temp),log(density), 1)
m = Const(1)
k = Const(2)

YBL = x.^m.*exp(k)

[Nicholas Green]

Cheers! Thanks!

[Star Strider]

My pleasure!

[Alessandro Saffioti]

@Star Strider. Great answer, can you explain how do you get this "x.^m.*exp(k)"?

[Star Strider]

Thank you.

A linear regression of the logarithms of the data produces:

log(y) = m*log(x) + k

so taking exp (the anti-logarithm) of both sides of the equation produces:

y = x^m * exp(k)

Note that taking the logarithms of the data distorts the errors (creating multiplicative errors from additive errors), so a nonlinear parameter estimation is more accurate.

Voting for my Answer would be appreciated.

[Alessandro Saffioti]

@Star Strider. If i'm using a semilogy plot how do i make the data linear? The result should be a natural logarithm

% code
r=cumsum([13 10 10 10 10 10])/1000;
temp=[30.0315 32.7117 31.4136 31.1503 29.5345 29.5345 29.3529];
y=polyval(polyfit(r,temp,2),r);
y2=polyval(polyfit(r,temp,1),r);
subplot(1,2,1)
plot(r,temp,'o:k',r,y)
subplot(1,2,2)
semilogy(r,temp, 'o:r',r,y2, 'k')

[Star Strider]

If your plot is:

log(y) = m*x + b

then:

y = exp(m*x + b)

or equivalently:

y = exp(b) * exp(m*x)

These are just basic log and exp identities.

I will let you adapt that to what you want to do.

[Marisabel Gonzalez]

@Star Strider how would you account for the errors using this method?

[Star Strider]

Good point.

I would now account for them by doing a direct parameter estimation of a nonlinear regression rather than a log-transformed linear regression.

Specifically:

temp = [100,200,400,600,800,1000,1200,1400,1600]';
density = [3.5,1.7,0.85,0.6,0.45,0.35,0.3,0.25,0.2]';

fcn = @(b,x) exp(b(1).*x).*exp(b(2)) + b(3);
[B,rsdnorm] = fminsearch(@(b) norm(density - fcn(b,temp)), [-0.01; max(density); min(density)]);
fprintf(1, 'Slope \t\t=%10.5f\nIntercept  \t=%10.5f\nOffset \t\t=%10.5f\n', B)

tv = linspace(min(temp), max(temp));

figure
plot(temp, density, 'p')
hold on
plot(tv, fcn(B,tv), '-')
grid

The errors are then the residuals. The norm of the residuals here is 0.342.

[P_L]

Hi Star Strider,

I have modified your code to fit my data.

I keep getting this error:

Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.594806 

Please could you help me fix this problem? Many Thanks

log_N = [-2.50000000000000,-2,-1.50000000000000,-1,-0.500000000000000,0,0.500000000000000]';
Log_R = [1.09861228866811,1.09861228866811,1.79175946922806,2.56494935746154,3.61091791264422,4.23410650459726,4.65396035015752]';
fcn = @(b,x) exp(b(1).*x).*exp(b(2)) + b(3);
[B,rsdnorm] = fminsearch(@(b) norm(Log_R - fcn(b,log_N)), [-0.01; max(Log_R); min(Log_R)]);
fprintf(1, 'Slope \t\t=%10.5f\nIntercept  \t=%10.5f\nOffset \t\t=%10.5f\n', B)
tv = linspace(min(log_N), max(log_N));
figure
plot(log_N, Log_R, 'p')
hold on
plot(tv, fcn(B,tv), '-')
grid

[Star Strider]

@P_L —

My pleasure.

Use an options structure:

opts = optimset('MaxFunEvals',5000);
[B,rsdnorm] = fminsearch(@(b) norm(Log_R - fcn(b,log_N)), [-0.01; max(Log_R); min(Log_R)], opts);

The rest of your code is unchanged.

[P_L]

Many many thanks :)

[Star Strider]

As always, my pleasure.

[Pepe]

Dear @Star Strider,

if T is in K and density is in kg/m^3, what would you say is the unit of the slope of the fitted line?

Is it (kg/m^3)/K or does it have log in it?

Thanks.

[Ivan Pixabaj]

Hey, I just wanted to ask that if you can modify the code to give estimated values for the density?

[Les Beckham]

Do it just as Star Strider did it above but with a different temperature

temp = [100,200,400,600,800,1000,1200,1400,1600]

temp = 1×9

100 200 400 600 800 1000 1200 1400 1600

density = [3.5,1.7,0.85,0.6,0.45,0.35,0.3,0.25,0.2]

density = 1×9

3.5000 1.7000 0.8500 0.6000 0.4500 0.3500 0.3000 0.2500 0.2000

%%Plot
loglog(temp, density, '*')
xlabel('Temperature T,K')
ylabel('Density D, kg/m^3')
title('Temperature vs. Density')
grid on

%%Bestfit Line
hold on
Const = polyfit(log(temp),log(density), 1)

Const = 1×2

-1.0023 5.8658

m = Const(1)

m = -1.0023

k = Const(2)

k = 5.8658

YBL = temp.^m.*exp(k)

YBL = 1×9

3.4901 1.7422 0.8697 0.5793 0.4342 0.3471 0.2892 0.2478 0.2167

loglog(temp,YBL)

densityAt500deg = 500.^m.*exp(k) % <<< Find a density estimate at a new temperature

densityAt500deg = 0.6954

plot(500, densityAt500deg, 'go') % <<< Plot the new estimate as a green circle

You could also use a vector of new temperatures instead of a scalar. Be careful of using points that aren't spanned by the range of the original fit. Extrapolation can be dangerous (give you misleading results).

Or adapt the non-linear regression approach instead.

[Ivan Pixabaj]

I get that but i was more asking along the lines of getting densitity values outside my density range, x range. For example let the extrapolation line not stop where my density stops but rather let it continue and estimate values outside the given range. Hopefully this clearly explains what i'm looking for, i appreciate the help. For instance let the line of extrapolation be continuous

[Star Strider]

@Ivan Pixabaj —

‘For example let the extrapolation line not stop where my density stops but rather let it continue and estimate values outside the given range.’

Do not ever extrapolate beyond the region of fit.

If you have a model for a process and you are reasonably certain that the available data describes the process and that the parameters estimated from it are reasonably reliable, it might be possible to extend the description of the process beyond the available data, providing the confidence intervals are also included.

If you only have random data and are doing curve fitting when the curve does not describe the actual process that created the data, this does not apply. You have absolutely no assurance that whatever created the available data will behave outside the limits of the data the same way it did within the limits of the data.

A simple example —

The nonlilnear parameter estimation approach is much more robust than fitting the logarithmically-transformed data, im part because it does not convert additive errors to multiplicative errors, as the logarithmic transformation does. So I am using it here.

Using the linear parameter estimation on logarithmically-transformed data would produce a similar result as this, however with a less accurate fit to the data.

x = logspace(0, 2, 150);
y = sin(2*pi*x/100)+1.5;
y = x.^4 - 14*x.^2 + 50;

fcn = @(b,x) x.^b(1).*exp(b(2));
B = fminsearch(@(b)norm(y(60:80) - fcn(b,x(60:80))), rand(2,1))

B = 2×1

4.3467 -0.9497

figure
loglog(x(60:80), y(60:80),'.-')
hold on
plot(x(60:80), fcn(B,x(60:80)), '-r')
hold off
grid
xlabel('x')
ylabel('y')

figure
loglog(x, y,'.')
hold on
plot(x, fcn(B,x), '-r', 'LineWidth',2)
hold off
grid
xline(x([60 80]), '--k', 'Limit of Fit')
xlabel('x')
ylabel('y')

So don’t extrapolate!

.

[Les Beckham]

Excellent example @Star Strider. Thumbs up.

I guess I should have emphasized that point more than with the cursory sentences "Be careful of using points that aren't spanned by the range of the original fit. Extrapolation can be dangerous (give you misleading results)."

[Star Strider]

@Les Beckham — Thank you!

I thought the concept needed an illustration, so I provided one.

Answer 2

[Walter Roberson]

logx = log(temp);
logy = log(density);
Const = polyfit(logx, logy, 1)
hold on
plot(temp, exp(polyval(Const, logx)));