logistic模型及其matlab算法

Logistic曲线的三种参数估计方法作者QQ：2377389590Logistic曲线的参数估计1844或1845年，比利时数学家PierreFrançoisVerhulst提出了logistic方程，这是一个对S型曲线进行数学描述的模型。一百多年来，这个方程多次应用于一些特殊的领域建模与预测，例如单位面积内某种生物的数量、人口数量等社会经济指标、某种商品（例如手机）的普及率等。logistic方程定义如下：1etbtxca(1)其中，t通常表示时间变量，,abc和为模型的参数；当趋势比较完整时0a0,0bc。其曲线如图1所示：0ln(c/a)/b01/(2c)1/c(ln(c/a)/b,1/(2c))根据1图和(1)方程式得：当t时，()0xt；当t时，()1/xtc。为了研究Logistic曲线的增长特性，对(1)式求导一阶导数得：20()btbtdxabedtcae图1logistic方程的曲线示意图设xt（t=1,2,…,n）为观测样本，对于logistic方程的参数，常规的估计方法有三种。1.1.Yule算法：根据式(1)，有111(1)=11()(1)(1)(1)tttttbtbtbtbbtbbtxxxxxcaecaeaeccecaeecex(2)设11ttttxxzx、1be以及(1)bce，则式(2)变形为线性方程ttzx，利用普通最小二乘（OLS）方法可以得到这个方程参数的估计值，b和c的估计值也可以进一步得到。为得到a的估计值，将式(1)变形为：1ˆˆˆlnlntcabtxt=1.2.…,n(3)左右分别对t求和11(1)ˆˆˆlnln2nttnncnabx(4)因此，a的估计值为：111(1)ˆˆˆexpln2nttnnacbnx(5)1.2.Rhodes算法：根据式（1），有(1)1(1)1(1)bttbbbtbbtcaexcceceaeecex(6)设11ttzx、1ttsx、(1)bce以及be，则式（6）变形为ttzs。用普通最小二乘（OLS）方法得到这个方程参数的估计值，并进一步得到b和c的估计值，然后利用式（3）-式（5）的方法得到a的估计值。1.3.Nair算法：式（2）可以进一步写为：11(1)(1)bttttxxecxx(7)即：11(1)11ttbttxcxexx(8)因此，有1111111121112(1)1211211bbbtttttttttttttteeexcxxxxcxxxxxcxxxx(9)式（9）可以整理为：11111112(1)11bbbbttttecexxexxe(10)设111tttzxx、111tttsxx、2(1)1bbcere以及11bbee，则式（10）可写为ttzs。用普通最小二乘（OLS）方法得到这个方程参数的估计值，并进一步得到b和c的估计值，然后利用式（3）-式（5）的方法得到a的估计值。2.列题：中国1965-2011年CO2排放量（单位：亿吨）如表1所示：表1中国1965-2011年CO2排放量1965196619671968196919701971197219731974480.9522.0468.8469.5573.8737.8869.8933.7977.2997.719751976197719781979198019811982198319841120.31176.11284.81422.11462.11499.71473.11539.21637.01771.019851986198719881989199019911992199319941886.51994.62145.72292.02396.82387.02484.42580.82750.22915.719951996199719981999200020012002200320043163.83231.93319.53319.63484.03550.63613.93833.14471.25283.020052006200720082009201020115803.26415.56797.97033.57636.38209.88979.1用logistic方程模拟我国CO2排放量的变化趋势，分别用三种方法估计方程参数，并分别计算三种方法的MAPE及未来五年CO2排放量的预测结果。2.1.Yule算法：clear;clc;%Yule算法：%author:朱伟杰%date:2018-1-24X=[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)-1fort=1:nZ(t)=(X(t+1)-X(t))/X(t+1)endX1=[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:2016YY=1./(c+a*exp(b*([XX-1965])))plot(XX,YY,'r-o')holdonplot(XX(1:length(X)),X,'g-^')legend('预测值','实际值')xlabel('年份');ylabel('CO_{2}排放量');title('CO_{2}预测值和实际值曲线图(Yule法)')set(gca,'XTick',[1965:2:2017])gridonformatshort;forecast=YY(end-4:end)%CO2排放量的预测结果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值a,b,c2.2.Rhodes算法：clear;clc;%Rhodes算法%author:朱伟杰%date:2018-1-24X=[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)-1fort=1:nZ(t)=1/X(t+1)S(t)=1/X(t)endX1=[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:2016YY=1./(c+a*exp(b*([XX-1965])))plot(XX,YY,'r-o')holdonplot(XX(1:length(X)),X,'k-^')set(gca,'XTick',[1965:2:2017])legend('预测值','实际值')xlabel('年份');ylabel('CO_{2}排放量');title('CO_{2}预测值和实际值曲线图(Rhodes法)')gridonformatshort;forecast=YY(end-4:end)%CO2排放量的预测结果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值a,b,c2.3.Nair算法：clear;clc;%Nair算法%author:朱伟杰%date:2018-1-24X=[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)-1fort=1:nZ(t)=1/X(t)-1/X(t+1)S(t)=1/X(t)+1/X(t+1)endX1=[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:2016YY=1./(c+a*exp(b*([XX-1965])))plot(XX,YY,'r-o')holdonplot(XX(1:length(X)),X,'g-^')legend('预测值','实际值')xlabel('年份');ylabel('二氧化碳排放量');title('二氧化碳预测值和实际值曲线图(Nair法)')set(gca,'XTick',[1965:2:2017])gridonformatshort;forecast=YY(end-4:end)%CO2排放量的预测结果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相对差值a,b,c根据编程计算，三种方法参数估计的结果及MAPE如表1所示：表1三种方法的参数估计结果及误差分析结果abcMAPEYule算法0.0020-0.0580-0.0000245790.1141Rhodes算法0.0019-0.08020.000109600.1259Nair算法0.0023-0.06830.0000165640.1271三种方法对未来五年CO2排放量的预测结果（单位：百万吨）如表2所示：表2中国CO2排放量的预测结果20122013201420152016Yule算法9