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

目录

 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方程定义如下：

                                                               x_{t}= \frac{1}{c+ae^{bt}}{\color{Red} }                                                         (1)

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

根据1图和(1)方程式得：当n \to +\infty时，x(t)\to\0 ；当t\to+\infty 时， x(t)\to\frac{1}{c}。为了研究Logistic曲线的增长特性，对(1)式求导一阶导数得：

                                                                \frac{\mathrm{d} x}{\mathrm{d} x}=\frac{-abe^{bt}}{(c+ae^{bt})^2}>0

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

1.1、Yule算法：

                                                                \dpi{120} \frac{x_{t+1}-x_{t}}{x_{t+1}}\\ =1-\frac{x_{t}}{x_{t+1}}\\ =1-\frac{c+ae^{b(t+1)}}{c+ae^{bt}}\\ =1-\frac{(ae^{bt}+c-c)(1-e^b)}{c+abe^{bt}}\\ =(1-e^b)-c(1-e^b)x_{t}                     (2)

根据式（1）有\dpi{120} Z_{t}=\frac{x_{t+1}-x_{t}}{x_{t+1}},\gamma =1-e^b，以及\beta =-(1-e^b)则式（2）变形为线性方程，Z_{t}=\gamma +\beta x_{t}利用普通最小二乘（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)  编辑  收藏  举报  来源

导航