Estimate r and K of logistic curve

I am given data that can be fitted to a logistic curve.

load population_america
data = [3.9290, 5.3080 ,7.2400, 9.6380, 12.8660, 17.0690, 23.1920, 31.4430, 38.5580, 50.1560, 62.9480, 75.9960, 91.9720, 105.7110, 122.7750,
    131.6690, 150.6970, 179.3230, 190.1850, 206.5460, 212.7100]
t = [ 1790:10:1990];

From this I want to estimate r and K and plot the observed data with the estimated function. I can assume that y0 = data(1) and it is stated that I should do linear regression on:

where z_i = 1/y_i and then use backslash to find the parameters r and K.

This is my code:

load population_america

y0=data(1);

% Plot observed data
figure(1);
plot(t,data,'*');
title('Exponential data fit');
xlabel('time');
ylabel('Population size');

% Make data linear
lnData = log(data);
% Plot linear data
figure(2)
plot(t,lnData,'*');
title('Linear data');
xlabel('time');
ylabel('ln(population size)');

% Estimate parameters using MATLAB row reduction
tData = [ones(size(t)),t];
parameters = tData\lnData;
r = parameters(1);
K = parameters(2);

% Generate data using the exponential model with estimated parameters
data_fit = exp(-r*t)*(1/y0)+(1-exp(-r*t))./K;
% Plot resulting data
figure(1);
hold on
plot(t,data_fit,'k-','linewidth',2);
legend('Observed data','Estimated function');

However, I don't know if this is correct nor does the estimated function show up in figure(1). Can someone spot the mistake? Thank you!

[VBBV]

Try this

 % if true
    % code
  %end
load population_america
y0=data(1);
% Make data linear
lnData = log(data);
% Plot linear data
figure(2)
plot(t,lnData,'*');
title('Linear data');
xlabel('time');
ylabel('ln(population size)');
% Estimate parameters using MATLAB row reduction
tData = [ones(size(t)),t];
parameters = tData\lnData;
r = parameters(1);
K = parameters(2);
% Generate data using the exponential model with estimated parameters
data_fit = exp(-r*t)*(1/y0)+(1-exp(-r*t))./K;
% Plot resulting data
figure(1);
plot(t,data,'*',t,data_fit,'k-','linewidth',2);
title('exponential data fit');
xlabel('time');
ylabel('population size');
legend('Observed data','Estimated function');

[katara]

That still doesn't work for me. The estimated function does not appear.

[katara]

I noticed that the function is wrong, if it is 1/y, then data_fit should be:

data_fit = 1./(exp(-r*t)*(1/data)+(1-exp(-r*t))./K);

But this still doesn't work.

I also saw that you could define a function as below and I tried using a different function to incorporate y0:

function y = logistic_func(t,K,r,y0)
y = (K*y0*exp(r.*t))./(y0*(exp(r.*t)-1)+K);
end

But then when I want to open it as:

data_fit = logistic_func

I get the error message:

Not enough input arguments.

Error in logistic_func (line 2)

y = (K*y0*exp(r.*t))./(y0*(exp(r.*t)-1)+K);

Error in logistic_curve (line 31)

data_fit = logistic_func;

Is there another way of opening the function?

Answer 1

[Alan Stevens]

Since y is your data you first need to invert the data to get z. Then you need to take ln(z) if you want to see a roughly straight line when plotted against t. You could then do a linear least squares fit to ln(z) vs t. You can extract the slope and intercept of this best fit line and compare the resut with the ln(z) values. However, because the logistic equation is a more complicated function of time you can't really extract r and K from this. However, here's a program that does the linear fit against ln(z).

% Population data ; y = population numbers, t = years
y = [3.9290, 5.3080 ,7.2400, 9.6380, 12.8660, 17.0690, 23.1920, 31.4430, ...
        38.5580, 50.1560, 62.9480, 75.9960, 91.9720, 105.7110, 122.7750, ...
        131.6690, 150.6970, 179.3230, 190.1850, 206.5460, 212.7100];
t = 1790:10:1990;
% Times need to be measured from 1790, so:
t = t - 1790;
% z values are inverted form of y
z = 1./y;
% Take logarithm of z
lnz = log(z);
% Now you have roughly linear relationship between z and t
% so do a least squares linear fit of lnz vs t
% i.e. look for; lnzlinear = m*t + c
M = [sum(t) numel(t); sum(t.^2) sum(t)];
V = [sum(lnz); sum(lnz.*t)];
mc = M\V;
m = mc(1);
c = mc(2);
% Now you can plot linear fit against data, where data is in the form 
% lnz vs t, i.e. log(1/y) vs t.
plot(t,lnz,'o',t,m*t+c),grid
xlabel('time'),ylabel('ln(z)')
legend('ln(1/y)','linear fit')

To get a reasonable estimate of r and K you need a more advanced fitting routine.

[katara]

This is my new code using the lnz from your code. However, I still get that the estimated function is a flat line.

load population_america


% Plot observed data
figure(1);
plot(t,data,'*');
title('Exponential data fit');
xlabel('time');
ylabel('Population size');


z = 1./data;
lnz = log(z);

% Estimate parameters using MATLAB row reduction
tData = [ones(size(lnz)),lnz];
parameters = tData\lnz;
b = parameters(1);
m = parameters(2);

% Tau = 10 so estimated parameters are
r = (-log(m)/10);
K=(1-exp(-r*10))/b;
y0=data(1);

% Generate data using the exponential model with estimated parameters

data_fit = (K.*y0.*exp(r.*t))./ (y0.*(exp(r.*t)-1)+K);

% Plot resulting data
figure(1);
hold on
plot(t,data_fit,'k-','linewidth',2);
legend('Observed data','Estimated function');

[Alan Stevens]

Look at the size of your variable “parameters”.

[katara]

Yeah I have that problem no matter how I change my z and lnz. I don't know why it doesn't work..

[Alan Stevens]

Try running my routine.

[katara]

This is the graph I get when I run your code. I don't know if this is correct however we have not used your type of coding:

M = [sum(t) numel(t); sum(t.^2) sum(t)];
V = [sum(lnz); sum(lnz.*t)];

And we are allowed to follow an example which uses for example, parameters = tData\lnData. That is why I wanted the code to work in that way and not use your code.

[Alan Stevens]

This is the curve that is produced by my routine:

But this is a linear fit to ln(1/y) and does not represent the logistic equation; it simply represents a straightforward exponential curve when turned back into y vs t. You need a more complicated fit to represent the logistic curve. Perhaps your statement " ... I should do linear regression ... " was not intended to mean that you wanted any sort of straight line fit to anything.

[Alan Stevens]

Having given this some more thought, here is a way of getting r and K from a linear regression approach. It assumes r*dt is small so that exp(-r*t) is expanded as 1 - r*dt between each data point. It's not a very good assumption here, but the following is the code and the resulting comparison graph:

% Population data ; y = population numbers, t = years
y = [3.9290, 5.3080 ,7.2400, 9.6380, 12.8660, 17.0690, 23.1920, 31.4430, ...
        38.5580, 50.1560, 62.9480, 75.9960, 91.9720, 105.7110, 122.7750, ...
        131.6690, 150.6970, 179.3230, 190.1850, 206.5460, 212.7100];
t = 1790:10:1990;
% Times need to be measured from 1790, so:
t = t - 1790;
% z values are inverted form of y
z = 1./y;

% z(i+1) = z(i)*exp(-r*t) + (1 - exp(-r*t))/K
% For small timeintervals: this becomes
% z(i+1) = z(i)*(1 - r*dt) + (r/K)*dt  or
% z(i+1) - z(i) = -r*z(i)*dt + (r/K)*dt  or
% dz(i) = [-z(i)*dt  dt]*[r; (r/K)] 
% Let M = [-z(i)*dt  dt]  so  M*[r; r/K] = dz
% Hence [r; r/K] = M\dz

dt = 10;                            % Constant timeintervals
dz = (z(2:end) - z(1:end-1))';      % Column vector of differences.
n = numel(dz);
% Linear regression
M = [-z(1:end-1)'  ones(n,1)]*dt;   % M is nx2 matrix
p = M\dz;                           % p = [r; r/K
r = p(1);
K = r/p(2);

% Plot comparison
tt = 0:t(end);
zfit = z(1).*exp(-r*tt) + (1 - exp(-r*tt))/K;
yfit = 1./zfit;

plot(t, y, 'o', tt, yfit), grid
xlabel('time'),ylabel('y')
legend('data','fit')
text(20,170,['r = ',num2str(r), '  K = ', num2str(K)])
% Not a good fit, because r*dt is not very small!