logistic模型及其matlab算法

1、Logistic曲線的三種參數(shù)估計(jì)方法作者QQ：2377389590Logistic曲線的參數(shù)估計(jì)1844或1845年，比利時(shí)數(shù)學(xué)家PierreFrangoisVerhulst提出了logistic方程，這是一個(gè)對(duì)S型曲線進(jìn)行數(shù)學(xué)描述的模型。一百多年來，這個(gè)方程多次應(yīng)用于一些特殊的領(lǐng)域建模與預(yù)測(cè)，例如單位面積內(nèi)某種生物的數(shù)量、人口數(shù)量等社會(huì)經(jīng)濟(jì)指標(biāo)、某種商品(例如手機(jī))的普及率等。logistic方程定義如下：1x(1)tc+aebt其中，t通常表示時(shí)間變量，a,b和c為模型的參數(shù)；當(dāng)趨勢(shì)比較完整時(shí)a0b0,c,0o其曲線如圖1所示：根據(jù)1圖和(1)方程式得：當(dāng)t-時(shí)，x(t)0；當(dāng)t+時(shí)，x

2、1/co為了研究Logistic曲線的增長(zhǎng)特性，對(duì)(1)式求導(dǎo)一階導(dǎo)數(shù)得：dx-abebt小,0dt(c+aebt)2圖1logistic方程的曲線示意圖設(shè)xt(t=1,2,，n)為觀測(cè)樣本，對(duì)于logistic方程的參數(shù)，常規(guī)的估計(jì)方法有三種。Yule算法：根據(jù)式(1)，有(7)t1xt1xt1caeb(t1)(2)caebt(aebt+c-c)(1-eb)caebt，(1-eb)-c(1-eb)xtx-xc設(shè)z二占F、丫，1-eb以及，-c（1-eb），則式變形為線性方程z，Y+x,txttt1利用普通最小二乘（OLS）方法可以得到這個(gè)方程參數(shù)的估計(jì)值，b和c的估計(jì)值也可以進(jìn)一步得到。為得

3、到a的估計(jì)值，將式（1）變形為：ln-c，ina+bt(二丿(3)(7)(7)左右分別對(duì)t求和1)-c，nIna+t=11xt丿a二exp廠1in1)-cn一t,11xt丿因此，a的估計(jì)值為：n(n1)b2n(n+1)b2(4)(5)(7)(7)Rhodes算法根據(jù)式（1），有1，c+aeb(t1)xt1，c-ceb+ceb+aeb(t+1)(6)eb，c(1-eb)+一x1設(shè)z=一txt11s，、丫，txttc（1-eb）以及，eb，則式（6）變形為z，丫+s。用普tt通最小二乘（OLS）方法得到這個(gè)方程參數(shù)的估計(jì)值，并進(jìn)一步得到b和c的估計(jì)值，然后利用式（3）-式（5）的方法得到a的估計(jì)值

4、。Nair算法：式（2）可以進(jìn)一步寫為：i+4_二(1-eb)(1-cx)xtt1(8)即：1x(1-cx)=t1t HYPERLINK l bookmark361一ebx一xt1t因此，有1+eb22x(1-cx)，+4t1xxt1tx2cxx+xt+1t+1tt(9)t+1t+1t+1式（9）可以整理為：111eb2c(1eb)1+eb(10)t+11111設(shè)z=、s，+txxtxxtt+1tt+1xxtt+12c(1eb)1ebr，以及卩，則式(10)可寫1+eb1+eb(8)(8)為z，丫+卩s。用普通最小二乘（OLS）方法得到這個(gè)方程參數(shù)的估計(jì)值，并進(jìn)一步得到ttb和c的估計(jì)值，然后

5、利用式（3）-式（5）的方法得到a的估計(jì)值。2.列題：中國(guó)1965-2011年CO2排放量（單位：億噸）如表1所示:表1中國(guó)1965-2011年CO2排放量1965196619671968196919701971197219731974480.9522.0468.8469.5573.8737.8869.8933.7977.2997.719751976197719781979198019811982198319841120.31176.11284.81422.11462.11499.71473.11539.21637.01771.019851986198719881989199019911992

6、199319941886.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方程模擬我國(guó)CO2排放量的變化趨勢(shì)，分別用三種方法估計(jì)方程參數(shù)，并分別計(jì)算三種方法的MAPE及未來五年CO

7、2排放量的預(yù)測(cè)結(jié)果。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,.3550.6,3613.9,3833.1,4471.2,5283,.5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1n=length(X)-1fort=1:nZ(t)=(X(t+1)-X(t)/X(t+1)endX1=ones(46,1)X(1:n)Y=ZB,Bint,r,rint,stats=regress(

8、Y,Xl)%最小二乘(0LS)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(預(yù)測(cè)值,實(shí)際值)xlabel(年份);ylabel(C0_2排放量);title(C0_2預(yù)測(cè)值和實(shí)際值曲線圖(Yule法)set(gca,XTick,1965:2:2017)gridonfor

9、matshort;forecast=YY(end-4:end)%C02排放量的預(yù)測(cè)結(jié)果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相對(duì)差值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,.3550.6,3613.9,3833.1,4471.2,5283,.n=length(X)-1fort=1:nZ(t)=1/X(t+1)S(t)=1/X(t)endX1=ones(46,

10、1)S(1:n)Y=ZB,Bint,r,rint,stats=regress(Y,Xl)%最小二乘(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(預(yù)測(cè)值,實(shí)際值)xlabel(年份);ylabel(

11、CO_2排放量);title(C0_2預(yù)測(cè)值和實(shí)際值曲線圖(Rhodes法)gridonformatshort;forecast=YY(end-4:end)%C02排放量的預(yù)測(cè)結(jié)果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相對(duì)差值a,b,c匚丐預(yù)測(cè)也和實(shí)麻由曲線圖很hod歸袪jmT1辨13711旳519771W11W1*9E忖怡11卿?S的陽2WJ?W5旳OS20n卻T左欄.8Nair算法：clear;clc;%Nair算法%author:朱偉杰%date:2018-1-24X=480.9,522,468.8,469.5,573.8,737.8,869

12、.8,933.7,977.2,.TOC o 1-5 h z.3550.6.3613.9.3833.1.4471.2.5283.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=ZB,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+

13、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(預(yù)測(cè)值,實(shí)際值)xlabel(年份);ylabel(二氧化碳排放量);title(二氧化碳預(yù)測(cè)值和實(shí)際值曲線圖(Nair法)set(gca,XTick,1965:2:2017)gridonformatshort;forecast=YY(end-4:end)%C02排放量的預(yù)測(cè)結(jié)果MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X)%平均相對(duì)差值a,b,c1WSS1!M|IMSIWrlIMiI-IJ.S-JJirIJUUb?Ett-d5WIWiWJ1J/J111TJ21舊匕IJ-J7TJLSh：I髭11ZD!-2Ph*一軌化軌預(yù)敵恒和實(shí)際恒曲繪圖小謝惟）1BKD5ZW-預(yù)瀏悄實(shí)際勺I根據(jù)編程計(jì)算，三種方法參數(shù)估計(jì)的結(jié)果及MAPE如表1所示：表1三種方法的參數(shù)估計(jì)結(jié)果及誤差分析結(jié)果abcMAPEYule算法0.0020-0.0580-0.0000245790.1141Rhodes算法0.0019-0.0802