Estimate r and K of logistic curve

11 次查看（过去 30 天）
katara
katara 2020-9-11
评论： Alan Stevens 2020-9-12
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 个评论
显示 1更早的评论
katara
katara 2020-9-11
That still doesn't work for me. The estimated function does not appear.
katara
katara 2020-9-11
编辑：katara 2020-9-11
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2020-9-11
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 个评论
显示 5更早的评论
Alan Stevens
Alan Stevens 2020-9-11
编辑：Alan Stevens 2020-9-11
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 2020-9-12
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!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by