田海麗

(鄭州工商學(xué)院 機(jī)械與電信工程學(xué)院，河南 鄭州451400)

基于MATLAB的二孩政策下人口增長模型研究
——以鄭州市人口為例

田海麗

(鄭州工商學(xué)院 機(jī)械與電信工程學(xué)院，河南 鄭州451400)

綜合考慮影響二孩的各種因素，以鄭州市歷年人口數(shù)據(jù)為基礎(chǔ)，給出二孩政策下人口增長模型函數(shù)——改進(jìn)的Logistic模型，并采用MATLAB進(jìn)行參數(shù)預(yù)測(cè)和函數(shù)建模；在建模過程中著重分析影響MATLAB建模性能優(yōu)劣的兩個(gè)關(guān)鍵因素：初值的設(shè)定和曲線擬合的實(shí)現(xiàn)方式；用MATLAB仿真軟件，嘗試多種Logistic曲線擬合的方法，找出性能良好的Logistic函數(shù)實(shí)現(xiàn)方法，最終確定二孩政策下的鄭州人口增長模型；對(duì)多種方法設(shè)置不同的初始值來檢驗(yàn)擬合方法對(duì)初始值的寬容度，最終找出對(duì)初始值寬容度較高的兩種MATLAB，實(shí)現(xiàn)Logistic曲線擬合的方法：nlinfit非線性擬合和lsqcurvefit最小二乘擬合.

二孩人口增長模型；MATLAB；Logic回歸模型；初始值設(shè)置；曲線擬合方法

0 引言

全面二孩政策放開以后，人口增加會(huì)對(duì)社會(huì)各方面帶來影響，人口增長模型研究迫在眉睫，根據(jù)鄭州市統(tǒng)計(jì)局公布的歷年人口數(shù)，選取1980—2015年鄭州市人口數(shù)據(jù)(表1)[1]進(jìn)行人口增長建模，并用MATLAB 2016a軟件仿真實(shí)現(xiàn).得到性能較好的人口增長模型，較為準(zhǔn)確地預(yù)測(cè)人口增長，進(jìn)而為社會(huì)的各項(xiàng)活動(dòng)順利開展提供人口數(shù)參考.

表1 鄭州歷年人口數(shù)Tab.1 The population of Zhengzhou in these years

1 相關(guān)理論和方法

1.1 Logic回歸模型[2]

荷蘭生物數(shù)學(xué)家提出Logistic人口增長模型為

(1)

其中，xm為自然資源和環(huán)境條件所能容許的最大人口，r為凈相對(duì)增長率(出生率減去死亡率)，x(t)為第t年的人口數(shù)，x0為初始人口數(shù).

考慮到影響二孩生育的因素，如夫妻年齡、夫妻職業(yè)、家庭收入、生育意愿、一孩性別、所在地區(qū)等，對(duì)Logistic人口增長模型進(jìn)行改進(jìn),得到

(2)

模型的關(guān)鍵是xm和r值的確定.

1.2 非線性回歸的MATLAB實(shí)現(xiàn)

1.2.1 MATLAB提供的函數(shù)實(shí)現(xiàn)非線性回歸

1)nlinfit函數(shù)[3],用于計(jì)算回歸系數(shù)，可以得出擬合的回歸系數(shù)和預(yù)測(cè)值.

2)lsqcurvefit函數(shù)[3],實(shí)現(xiàn)非線性最小二乘擬合[3]，命令格式為：[x,resnorm,res]=lsqcurvefit(@f,x0,t,y,lb,ub),f可以是M-函數(shù)或inline函數(shù)，x0=參數(shù)初值；t=自變量向量； y=函數(shù)值向量；lb,ub=參數(shù)的上下界；輸出 x=參數(shù)估計(jì)值；resnom=(誤差平方和)；res=r(誤差向量).

3)lsqnonlin函數(shù)[4],實(shí)現(xiàn)非線性最小二乘擬合,命令格式為x=lsqnonlin(fun，x0，lb,ub)：lb與ub為變量x的下界及上界.lsqnonlin函數(shù)采用的是迭代法，由于程序的局限性，不可能搜索無窮大的區(qū)間，因此初始值的選擇很重要.如果最優(yōu)解離所給初始值比較近，迭代出該最優(yōu)解的概率就很高.lsqnonlin指令適合求解任意形式的非線性擬合，需要先定義擬合函數(shù)，任取初始向量，調(diào)用lsqnonlin命令即可求得最優(yōu)值.

1.2.2 MATLAB自定義非線性回歸函數(shù)

MATLAB可以用3種方式來自定義函數(shù)：①命令inline;②用 “.m”文件；③用匿名函數(shù)定義要擬合的函數(shù).

1.2.3 非線性回歸曲線擬合相關(guān)方法和關(guān)鍵

1)初值的設(shè)定.非線性回歸通常的求解方法具有局部收斂性，需要給出合適的初值，從生態(tài)學(xué)角度考慮，初始種群大小對(duì)于種群增長非常重要，初值準(zhǔn)確與否是整個(gè)求解工作成敗的關(guān)鍵[2,5-6]，因此需要找出對(duì)初值容忍度較大的實(shí)現(xiàn)方法.

2)曲線擬合實(shí)現(xiàn)方法.①非線性擬合模型線性化,常用的有取變量的倒數(shù)、對(duì)數(shù)，關(guān)系式的微分-積分線性化變化，分段線性化，將非線性擬合問題轉(zhuǎn)化為線性擬合問題[5]，線性化擬合函數(shù)常用的有polyfit函數(shù)和regress線性回歸函數(shù).②直接非線性化擬合,常用的有nlinfit擬合、lsqcurvefit函數(shù)擬合、lsqnonlin函數(shù)擬合.

1.3 曲線擬合度的評(píng)判準(zhǔn)則

決定系數(shù)R2，也稱為擬合優(yōu)度，R2越接近于1越好，越接近1表示擬合度越好，R2的公式為

2 MATLAB曲線擬合實(shí)現(xiàn)

2.1 非線性直接擬合

1)取倒數(shù)，用nlinfit函數(shù)進(jìn)行非線性擬合

實(shí)現(xiàn)程序如下，本次程序完整可運(yùn)行，后面只給出程序核心部分.

>> t=0:1:35;

x=[450.6 458.8 467.6 474.0 480.0 485.3 491.5 500.5 510.7 521.4 557.8 565.1 570.9 581.2 588.6 600.3 607.6 614.8 622.7 631.6 665.9 677.0 687.7 697.7 708.2 716.0 724.3 735.6 743.6 752.1 866.1 885.7 903.1 919.1 937.8 956.9];

y=1./x; %取倒數(shù)

fy=inline(′a(1)+(1/450.6-a(1))*exp(-1.62*a(2)*t)′,′a′,′t′) %定義擬合函數(shù)

a=nlinfit(t,y,fy,[200,0.7]) %進(jìn)行非線性擬合

t1=min(t):0.1:max(t);

y1=a(1)+(1/450.6-a(1))*exp(-1.62*a(2)*t1);

y2=1./y1; %還原

xm=1/a(1)

r=a(2)

y3=fy(a,t);

y4=1./y3;

R2=1-sum((y-y3).^2)/sum((y-mean(y)).^2) %決定系數(shù)

plot(t,x,′*′,t,y4,′o′,t1,y2)

>> title(′Logistic擬合曲線′)

>> legend(′Data′,′Logistic Fit′)

此時(shí)，決定系數(shù)為R2=0.989 1，接近于1，擬合優(yōu)度較高，人口模型函數(shù)為

(3)

這種方法在確定a=nlinfit(t,y,fy,[-1,1])時(shí)，對(duì)初始值取值的寬容度較大，具體數(shù)據(jù)見圖1，擬合結(jié)果見表2.

圖1 nlinfit函數(shù)進(jìn)行非線性擬合Fig.1 The performance of nlinfit nonlinear fitting

表2 不同初始值擬合結(jié)果Tab.2 Effect of different initial values

2)lsqcurvefit最小二乘擬合實(shí)現(xiàn)

主要程序如下：

x0=[500,0.02];

[x,norm,res]=lsqcurvefit(@mylogistic,x0,t,y);

y1=mylogistic(x,t);

R2=0.9666

而且可以驗(yàn)證這種方法對(duì)初值選擇有比較大的容忍度，且容易實(shí)現(xiàn)，擬合系數(shù)R2接近于1，具體數(shù)據(jù)見圖2，擬合結(jié)果見表3.

3)lsqnonlin求最優(yōu)解擬合

最小二乘法(least square method)，通過最小誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配，可用于曲線擬合.在MATLAB 5.x中，用函數(shù)leastsq[8]，在MATLAB 2016a中使用lsqnonlin.

function f=cf1(a)

t=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35];

y=[450.6,458.8,467.6,474.0,480.0,485.3,491.5,500.5,510.7,521.4,557.8,565.1,570.9,581.2,588.6,600.3,607.6,614.8,622.7,631.6,665.9,677.0,687.7,697.7,708.2,716.0,724.3,735.6,743.6,752.1,866.1,885.7,903.1,919.1,937.8,956.9];

f=a(1)/(1+(a(1)/450.6-1))*exp(-1.62*a(2)*t);

end

圖2 lsqcurvefit最小二乘擬合Fig.2 The performance of lsqcurvefit nonlinear fitting

表3 不同初始值擬合結(jié)果Tab.3 Effect of different initial values

以文件名cf1.m存盤(文件名與函數(shù)名應(yīng)一致)，然后在MATLAB的工作窗口輸入語句： a0=[1,1];a=lsqnonlin (@cf1,a0);成功取得最優(yōu)解為：a=[331.030 9,5.265 7]，R2= -3.055 1，擬合結(jié)果較差.進(jìn)行其他多次嘗試,都沒有好的擬合結(jié)果(表4).這種方法最優(yōu)解的收斂與初值選擇很有關(guān)系，初值選擇不好，擬合出的參數(shù)結(jié)果較差，試解很多次，仍然找不出適合的參數(shù)值.此命令在進(jìn)行二孩人口建模時(shí)不適用.

表4 不同初始值擬合結(jié)果

2.2 非線性曲線線性化擬合

1)求導(dǎo)后用polyfit進(jìn)行多項(xiàng)式線性擬合

主要程序如下：

>> x0=450.6;

>>y=diff([x0 x])./(1*x);%一階求導(dǎo)

>> A=polyfit(x,y,1);%多項(xiàng)式線性曲線擬合

>> polyval(A,t)%多項(xiàng)式在x處的擬合值y

>> r=A(:,2)；s=-A(:,1)；xm=r/s;

此時(shí)，xm = 95.646 8，r =-0.003 5，此時(shí)R2 =0.911 8.人口模型參數(shù)為

(4)

擬合圖與原始圖對(duì)比曲線見圖3.

2)線性化后regress回歸

將dx/dt=rx(1-x/xm)改為：(dx/dt)/x=r-sx,其中s=r/xm這樣對(duì)r,s的估計(jì)是線性的，左邊可用數(shù)值微分計(jì)算.主要程序如下：

>> dx=diff(x);dxx=dx./x(1:35); dxx=dxx′; tt=x(1:35); z=[ones(35,1) tt′];

>> b=regress(dxx,z);z(1)=-b(1)/b(2);z(2)=b(1);

>> y1=mylogistic(z,t);

由圖4可知，擬合結(jié)果較差，不宜采用此方法.

圖3 求導(dǎo)后多項(xiàng)式polyfit擬合Fig.3 The performance of polyfit linear fitting after derivation

圖4 線性化后regress回歸Fig.4 The performance of regress after linearization

3 結(jié)論

在用MATLAB做非線性回歸時(shí)，實(shí)現(xiàn)的方法有很多種，而且不同的初始值和不同的實(shí)現(xiàn)方法會(huì)影響到擬合函數(shù)是否可有效利用，需要進(jìn)行不斷地嘗試.可得出以下結(jié)論：

1) 曲線非線性化直接擬合結(jié)果優(yōu)于曲線線性化擬合結(jié)果.

2) 用nlinfit函數(shù)進(jìn)行非線性Logsitic改進(jìn)模型擬合時(shí)，在確定系數(shù)時(shí)，對(duì)初值的容忍度較大，擬合系數(shù)R2=0.989 1，擬合性能優(yōu)于lsqcurvefit最小二乘擬合，是實(shí)現(xiàn)方法中性能最好的.

3) lsqcurvefit最小二乘擬合進(jìn)行Logsitic改進(jìn)模型擬合時(shí)，對(duì)初值設(shè)置的寬容度較大，比較容易實(shí)現(xiàn)，擬合系數(shù)R2=0.966 6，雖然擬合系數(shù)稍劣于nlinfit擬合，但是R2接近于1，由此方法確定的函數(shù)模型可以接受.

4) lsqnonlin函數(shù)在其他函數(shù)建模中表現(xiàn)出較好的性能，但是在二孩人口建模時(shí)找的最優(yōu)解擬合性能很差，無法采用.

[1] 鄭州市統(tǒng)計(jì)局.歷年人口數(shù)[EB/OL]. 2017-01-01[2017-02-09].http://tjj.zhengzhou.gov.cn/ndsj/216707.jhtml.

[2] 楊益民,付必謙.關(guān)于Logistic 增長模型參數(shù)估計(jì)方法的再探討[J].統(tǒng)計(jì)與決策,2015(13):28-29.

[3] 馬莉.MATLAB數(shù)學(xué)實(shí)驗(yàn)與建模[M].北京：清華大學(xué)出版社，2010:148-149.

[4] 何正風(fēng).MATLAB在數(shù)學(xué)方面的應(yīng)用[M]. 北京：清華大學(xué)出版社，2012:225-226.

[5] 胡亮. 非線性擬合的初值問題[J]. 吉首大學(xué)學(xué)報(bào) (自然科學(xué)版), 2003, 24(1): 37-39.

[6] 殷祚云. Logistic 曲線擬合方法研究[J]. 數(shù)理統(tǒng)計(jì)與管理, 2002, 21(1): 41-46.

[7] 謝蘭,高東紅. 非線性回歸方法的應(yīng)用與比較[J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí), 2009, 39(10): 117-121.

[8] 袁新生, 龍門. 非線性曲線擬合的三種軟件解法之比較[J]. 徐州工程學(xué)院學(xué)報(bào), 2005, 20(3): 74-76.

ResearchedonModelofPopulationGrowthUnderTwo-childPolicyBasedonMATLAB：ACaseStudyofZhengzhouCity

TIAN Haili
(CollegeofMechanicalandTelecommunicationsEngineering,ZhengzhouTechnologyandBusinessUniversity，Zhengzhou451400,China)

Taking factors on bearing second-child into account, based on data of Zhengzhou City in recent years, an population growth model, namely, an improved logistic population model, under two-child policy was presented. Then the parameters were estimated and functions were modeled based on MATLAB 2016a. It focused on two chief factors affecting the performance of curve fitting method with MATLAB: setting of initial parameter value and curve fitting methods. After repeated attempts in variety of implementations with MATLAB，the better method of curve fitting was founded, and the population growth model of Zhengzhou under two-child policy was established finally. Moreover, tested the different nonlinear fitting methods’ tolerance of initial value, and found the best two methods, namely, the nlinfit nonlinear fitting method and the lsqcurvefit least square fitting method.

the population growth model under the two-child policy; MATLAB; logistic regress model; setting of initial parameter value; curve fitting method

2017-04-13

河南省社科聯(lián)、河南省經(jīng)團(tuán)聯(lián)調(diào)研課題“基于二胎政策下人口增長模型對(duì)鄭州市中小學(xué)教育發(fā)展的影響”(SKL-2016-697)；2015年河南省教育技術(shù)裝備和實(shí)踐教育研究項(xiàng)目“高職移動(dòng)通信技術(shù)專業(yè)公共實(shí)踐教學(xué)平臺(tái)建設(shè)研究”(GZS222)

田海麗(1978—)，女，河南開封人，鄭州工商學(xué)院機(jī)械與電信工程學(xué)院講師.

10.3969/j.issn.1007-0834.2017.02.008

TP391

：A

：1007-0834(2017)02-0030-06