Estimate r and K of logistic curve

13 visualizaciones (últimos 30 días)
katara
katara el 11 de Sept. de 2020
Comentada: Alan Stevens el 12 de Sept. de 2020
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!
  3 comentarios
katara
katara el 11 de Sept. de 2020
That still doesn't work for me. The estimated function does not appear.
katara
katara el 11 de Sept. de 2020
Editada: katara el 11 de Sept. de 2020
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?

Iniciar sesión para comentar.

Respuestas (1)

Alan Stevens
Alan Stevens el 11 de Sept. de 2020
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.
  7 comentarios
Alan Stevens
Alan Stevens el 11 de Sept. de 2020
Editada: Alan Stevens el 11 de Sept. de 2020
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
Alan Stevens el 12 de Sept. de 2020
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!

Iniciar sesión para comentar.

Community Treasure Hunt

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

Start Hunting!

Translated by