Logistic三种迭代算法公式推导及matlab代码实现_logistic方程

目录

 Logistic曲线的参数估计

1.1、Yule算法：

1.2、Rhodes算法：

1.3、Nair算法：

2、matlab算法实现

2.1.Yule算法:

2.2.Rhodes算法：

2.3.Nair算法：

点击链接加入群聊【matlab&&Python划水群】

---

---

         

 Logistic曲线的参数估计

作者qq:2377389590欢迎探讨。

闲着无聊整理了一下课程设计的作业，老师给了我们数学模型，让我们根据数学模型原理写代码。聪明的小编，在此夸夸自己。废话少说，原理如下：

       1844或1845年，比利时数学家Pierre François Verhulst提出了logistic方程，这是一个对S型曲线进行数学描述的模型。一百多年来，这个方程多次应用于一些特殊的领域建模与预测，例如单位面积内某种生物的数量、人口数量等社会经济指标、某种商品（例如手机）的普及率等。

logistic方程定义如下：

                                                                                                                        (1)

其中，t通常表示时间变量， 为模型的参数；当趋势比较完整时   a>0  b<0  c>0；其曲线如图1所示：

根据1图和(1)方程式得：当时， ；当 时， 。为了研究Logistic曲线的增长特性，对(1)式求导一阶导数得：

                                                                

设xt（t=1, 2, …, n）为观测样本，对于logistic方程的参数，常规的估计方法有三种。

1.1、Yule算法：

                                                                                     (2)

根据式（1）有,，以及则式（2）变形为线性方程，利用普通最小二乘（OLS）方法可以得到这个方程参数的估计值，b和c的估计值也可以进一步得到。

为得到a的估计值，将式（1）变形为：

               

1.2、Rhodes算法：

根据式（1），有

普通最小二乘（OLS）方法得到这个方程参数的估计值，并进一步得到b和c的估计值，然后利用式（3）-式（5）的方法得到a的估计值

1.3、Nair算法：

（2）式可以进一步写为：

2、matlab算法实现

中国1965-2011年CO2排放量（单位：亿吨）如表1所示：(备注：这是上课老师给的数据)

表1中国1965-2011年CO2排放量

用logistic方程模拟我国CO2排放量的变化趋势，分别用三种方法估计方程参数，并分别计算三种方法的MAPE及未来五年CO2排放量的预测结果。

2.1.Yule算法:

                       

clear;clc;
%Yule算法
%E-mail：2377389590@qq.com
%date:2018-1-24
X=[480.9,522,468.8,469.5,573.8,737.8,869.8,933.7,977.2,...
    997.7,1120.3,1176.1,1284.8,1422.1,1462.1,1499.7,...
    1473.1,1539.2,1637,1771,1886.5,1994.6,2145.7,2292,...
    2396.8,2387,2484.4,2580.8,2750.2,2915.7,3163.8,3231.9,...
    3319.5,3319.6,3484.,3550.6,3613.9,3833.1,4471.2,5283,...
    5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1]
n=length(X)-1
for t=1:n
    Z(t)=(X(t+1)-X(t))/X(t+1)
end
X1=[ones(46,1) X(1:n)']
Y=Z'
[B,Bint,r,rint,stats]=regress(Y,X1)% 最小二乘（OLS）
gamma=B(1,1)
beta=B(2,1)
b=log(1-gamma)
c=beta/(exp(b)-1)
 a=exp((sum(log(1./X(1:n)-c))-n*(n+1)*b/2)/n)
XX=1965:2016
YY=1./(c+a*exp(b*([XX-1965])))
plot(XX,YY,'r-o')
hold on
plot(XX(1:length(X)),X,'g-^')
legend('预测值','实际值')
xlabel('年份');ylabel('CO_{2}排放量');
title('CO_{2}预测值和实际值曲线图(Yule法)')
set(gca,'XTick',[1965:2:2017])
grid on
format short;
forecast=YY(end-4:end)%CO2排放量的预测结果
MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值
a,b,c

            

2.2.Rhodes算法：

clear;clc;
%Rhodes算法：
%E-mail：2377389590@qq.com
%date:2018-1-24
X=[480.9,522,468.8,469.5,573.8,737.8,869.8,933.7,977.2,...
    997.7,1120.3,1176.1,1284.8,1422.1,1462.1,1499.7,...
    1473.1,1539.2,1637,1771,1886.5,1994.6,2145.7,2292,...
    2396.8,2387,2484.4,2580.8,2750.2,2915.7,3163.8,3231.9,...
    3319.5,3319.6,3484.,3550.6,3613.9,3833.1,4471.2,5283,...
    5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1]
n=length(X)-1
for t=1:n
    Z(t)=1/X(t+1)
    S(t)=1/X(t)
end
X1=[ones(46,1) S(1:n)']
Y=Z'
[B,Bint,r,rint,stats]=regress(Y,X1)% 最小二乘（OLS）
gamma=B(1,1)
beta=B(2,1)
b=log(beta)
 c=gamma/(1-exp(b))
 a=exp((sum(log(1./X(1:n+1)-c))-(n+1)*(n+2)*b/2)/(n+1))
XX=1965:2016
YY=1./(c+a*exp(b*([XX-1965])))
plot(XX,YY,'r-o')
hold on
plot(XX(1:length(X)),X,'k-^')
set(gca,'XTick',[1965:2:2017])
legend('预测值','实际值')
xlabel('年份');ylabel('CO_{2}排放量');
title('CO_{2}预测值和实际值曲线图(Rhodes法)')
grid on
format short;
forecast=YY(end-4:end)%CO2排放量的预测结果
MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值
a,b,c

2.3.Nair算法：

clear;clc;
%Nair算法 
%E-mail：2377389590@qq.com
%date:2018-1-24
X=[480.9,522,468.8,469.5,573.8,737.8,869.8,933.7,977.2,...
    997.7,1120.3,1176.1,1284.8,1422.1,1462.1,1499.7,...
    1473.1,1539.2,1637,1771,1886.5,1994.6,2145.7,2292,...
    2396.8,2387,2484.4,2580.8,2750.2,2915.7,3163.8,3231.9,...
    3319.5,3319.6,3484.,3550.6,3613.9,3833.1,4471.2,5283,...
    5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1]
n=length(X)-1
for t=1:n
    Z(t)=1/X(t)-1/X(t+1)
    S(t)=1/X(t)+1/X(t+1)
end
X1=[ones(46,1) S(1:n)']
Y=Z'
[B,Bint,r,rint,stats]=regress(Y,X1)% 最小二乘（OLS）
gamma=B(1,1)
beta=B(2,1)
b=log((1-beta)/(1+beta))
c=gamma*(1+exp(b))/(2*(exp(b)-1))
a=exp((sum(log(1./X(1:n)-c))-n*(n+1)*b/2)/n)
XX=1965:2016
YY=1./(c+a*exp(b*([XX-1965])))
plot(XX,YY,'r-o')
hold on
plot(XX(1:length(X)),X,'g-^')
legend('预测值','实际值')
xlabel('年份');ylabel('二氧化碳排放量');
title('二氧化碳预测值和实际值曲线图(Nair法)')
set(gca,'XTick',[1965:2:2017])
grid on
format short;
forecast=YY(end-4:end)%CO2排放量的预测结果
MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值
a,b,c

根据编程计算，三种方法参数估计的结果及MAPE如表1所示：

三种方法的参数估计结果及误差分析结果

三种方法对未来五年CO2排放量的预测结果（单位：百万吨）如表2所示：

 中国CO2排放量的预测结果

如果想要最原始的文件和代码可以查看我上传的资源，源码下载地址

转载请备注

点击链接加入群聊【matlab&&Python划水群】