Logistic三种迭代算法公式推导及matlab代码实现
目录
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排放量
1965 | 1966 | 1967 | 1968 | 1969 | 1970 | 1971 | 1972 | 1973 | 1974 |
480.9 | 522.0 | 468.8 | 469.5 | 573.8 | 737.8 | 869.8 | 933.7 | 977.2 | 997.7 |
1975 | 1976 | 1977 | 1978 | 1979 | 1980 | 1981 | 1982 | 1983 | 1984 |
1120.3 | 1176.1 | 1284.8 | 1422.1 | 1462.1 | 1499.7 | 1473.1 | 1539.2 | 1637.0 | 1771.0 |
1985 | 1986 | 1987 | 1988 | 1989 | 1990 | 1991 | 1992 | 1993 | 1994 |
1886.5 | 1994.6 | 2145.7 | 2292.0 | 2396.8 | 2387.0 | 2484.4 | 2580.8 | 2750.2 | 2915.7 |
1995 | 1996 | 1997 | 1998 | 1999 | 2000 | 2001 | 2002 | 2003 | 2004 |
3163.8 | 3231.9 | 3319.5 | 3319.6 | 3484.0 | 3550.6 | 3613.9 | 3833.1 | 4471.2 | 5283.0 |
2005 | 2006 | 2007 | 2008 | 2009 | 2010 | 2011 |
|
|
|
5803.2 | 6415.5 | 6797.9 | 7033.5 | 7636.3 | 8209.8 | 8979.1 |
用logistic方程模拟我国CO2排放量的变化趋势,分别用三种方法估计方程参数,并分别计算三种方法的MAPE及未来五年CO2排放量的预测结果。
2.1.Yule算法:
clear;clc;
%Yule算法
%author:zhuweijie
%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算法:
%author:zhuweijie
%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算法
%author:zhuweijie
%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所示:
三种方法的参数估计结果及误差分析结果
a | b | c | MAPE | |
Yule算法 | 0.0020 | -0.0580 | -0.000024579 | 0.1141 |
Rhodes算法 | 0.0019 | -0.0802 | 0.00010960 | 0.1259 |
Nair算法 | 0.0023 | -0.0683 | 0.000016564 | 0.1271 |
三种方法对未来五年CO2排放量的预测结果(单位:百万吨)如表2所示:
中国CO2排放量的预测结果
2012 | 2013 | 2014 | 2015 | 2016 | |
Yule算法 | 9167 | 9847 | 10587 | 11396 | 12282 |
Rhodes算法 | 6476.7 | 6624.9 | 6767.8 | 6905.3 | 7037.3 |
Nair算法 | 9050 | 9588 | 10152 | 10742 | 11359 |
如果想要最原始的文件和代码可以查看我上传的资源,源码下载地址
点击链接加入群聊【matlab&&Python划水群】
posted on 2019-11-25 09:45 好玩的MATLAB 阅读(27) 评论(0) 编辑 收藏 举报 来源