馬爾薩斯數(shù)據(jù)擬合實(shí)驗(yàn)

1、馬爾薩斯數(shù)據(jù)擬合實(shí)驗(yàn)實(shí)驗(yàn)內(nèi)容實(shí)驗(yàn)一：數(shù)據(jù)擬合Malthus人口指數(shù)增長(zhǎng)模型中的參數(shù)1790-1980年間美國(guó)每隔10年的人口記錄如表5.3：表5.3年份1790180018101820183018401850人口（ X106）3.95.37.29.612.917.123.2年份1860187018801890190019101920人口（ X106）31.438.650.262.976.092.0106.5年份193019401950196019701980人口（ X10 6）123.2131.7150.7179.3204.0226.5用以上數(shù)據(jù)檢驗(yàn)馬爾薩斯（Malthus）人口指數(shù)增長(zhǎng)模型，

2、根據(jù)檢驗(yàn)結(jié)果進(jìn)一步討論馬爾薩斯人口模型 的改進(jìn)。試驗(yàn)二：經(jīng)濟(jì)增長(zhǎng)模型增加生產(chǎn)、發(fā)展經(jīng)濟(jì)所依靠的主要因素有增加投資、增加勞動(dòng)力以及技術(shù)革新等，在研究國(guó)名經(jīng)濟(jì)產(chǎn)值 與這些因素的數(shù)量關(guān)系時(shí)，由于技術(shù)水平不像資金、勞動(dòng)力那樣容易定量化，作為初步的模型，可認(rèn)為技 術(shù)水平不變，只討論不像資金、勞動(dòng)力之間的關(guān)系。在科學(xué)技術(shù)發(fā)展不快時(shí)，如資本主義經(jīng)濟(jì)發(fā)展的前期， 這種模型是有意義的。用Q,K,L分別表示產(chǎn)值、資金、勞動(dòng)力，要尋求的數(shù)量關(guān)系Q（K，L）。經(jīng)過(guò)簡(jiǎn)化假設(shè)與分析，在經(jīng)濟(jì) 學(xué)中，推導(dǎo)出一個(gè)著名的Cobb-Douglas生產(chǎn)函數(shù)Q（K, L ）= aK aLP,0 , P 1（*）式中a, P,a要由緊急

3、統(tǒng)計(jì)數(shù)據(jù)確定?，F(xiàn)有美國(guó)馬薩諸賽州1900-1926年上述三個(gè)經(jīng)濟(jì)指數(shù)的統(tǒng)計(jì)數(shù)據(jù)。如表5.5,試用數(shù)據(jù)擬合的方法，求出（*）式中的參數(shù)a, P,a。表5.5tQKLtQKL19001.051.041.0519142.013.241.6519011.181.061.0819152.003.241.6219021.291.161.1819162.092.611.8619031.301.221.2219171.964.101.9319041.301.271.1719182.204.361.9619051.421.371.3019192.124.771.9519061.501.441.3919202.1

4、64.751.9019071.521.531.4719212.084.541.5819081.461.571.3119222.244.541.6719091.602.051.4319232.564.581.8219101.692.511.5819242.344.581.6019111.812.631.5919252.454.581.6119121.932.741.6619262.584.541.6419131.952.821.68三、實(shí)驗(yàn)環(huán)境Windows操作系統(tǒng)；MATLAB 7.0.四、實(shí)驗(yàn)過(guò)程一、問(wèn)題：實(shí)驗(yàn)一二、問(wèn)題分析因?yàn)榧僭O(shè)人口增長(zhǎng)率為常數(shù)，且已知初始時(shí)刻的人口數(shù)，可根據(jù)t時(shí)刻人口變

5、化用差商代替微商。當(dāng) 菩T 0時(shí)，可列出關(guān)于t時(shí)刻人口變化的微分方程。dx=TX dt。根據(jù)該微分方程求出解，再用數(shù)據(jù)擬合模型球參數(shù)。x（0 ）二 xI 0三、假設(shè)1）假設(shè)人口增長(zhǎng)率為常數(shù)，記為r.2）記時(shí)刻t的人口為x（t）（即x（t）為模型的狀態(tài)變量）。3）假設(shè)初始時(shí)刻人口為x0。四、模型建立1）求解該微分方程得到x（t）的表達(dá)式為xD = xo *e，*t則可轉(zhuǎn)化為求關(guān)于x（t）在t時(shí)刻的測(cè)量數(shù)據(jù)來(lái)辨識(shí)參數(shù)對(duì)應(yīng)的數(shù)學(xué)模型球函數(shù)的最小二乘解。其中a=0,b=x0，則需求r的最小值點(diǎn)。2）運(yùn)用最小二函數(shù)擬合求出所需的r值。五、模型求解編寫 M 文件 curvefun.mfunction f=

6、curvefun(x,tdata)f=1e+05*39*exp(x*tdata);編寫(testl.m)tdata=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;x0=0.2;y=lsqcurvefit(curvefun,x0,tdata,cdata)即求出結(jié)果為：y = 0.0211即 r= 0.0211根據(jù)下述程序作出模型圖像：tdata=linspace(0,200,20);cdata=1e+05.*39

7、53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;plot(tdata,cdata)hold onf=1e + 05*39*exp(0.0211 *tdata);plot(tdata,f,r),title( 馬爾薩斯模型)模型圖像：六、模型評(píng)價(jià)可以看出使用一般曲線擬合，出來(lái)走勢(shì)基本一致，但是還有一定誤差，并且由圖可以看出該曲線接近一條開(kāi)口向上的二次曲線。因此可以嘗試使用二次多項(xiàng)式對(duì)馬爾薩斯模型進(jìn)行擬合。將會(huì)大大簡(jiǎn)化求解的過(guò)程。七、模型改進(jìn)即推廣則二次曲線擬合求解如下：編寫程序：tdat

8、a=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;p=polyfit(tdata,cdata,2)得到結(jié)果：p =1.0e+006 *0.0059-0.08265.4634編程作圖：tdata=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2

9、265;plot(tdata,cdata)hold onf=1.0e+006 * 5.4634.*tdata.A2-1.0e+006 *0.0826.*tdata+1.0e+006 * 0.0059；plot(tdata,f,r),title( 馬爾薩斯模型 1)則可根據(jù)上述程序作出二次曲線與點(diǎn)線圖如下：可見(jiàn)該二次曲線與馬爾薩斯模型基本上是完全重合的。因此兒二次曲線擬合是優(yōu)于一般的曲線擬合 的，對(duì)于馬爾薩斯模型來(lái)說(shuō)。、問(wèn)題：實(shí)驗(yàn)二二、問(wèn)題分析因?yàn)橐阎a(chǎn)函數(shù)為Q(K L，且已知初始時(shí)刻的相關(guān)Q、K、L，可根據(jù)上式進(jìn)行數(shù)據(jù)擬合，但是需將其轉(zhuǎn)化為2維數(shù)組計(jì)算即將L、K合成為一個(gè)2維數(shù)組data再?gòu)?/p>

10、data中將K、L取出， 則可利用一維擬合求出所要的解。三、假設(shè)Q,K,L分別表示產(chǎn)值、資金、勞動(dòng)力認(rèn)為技術(shù)水平不變，只討論不像資金、勞動(dòng)力之間的關(guān)系四、模型建立五、模型求解編程如下：編寫M文件curveu.mfunction q=curvev(x,data)K=data(1,:);L=data(2,:);q=x(1).*(K.八x(2).*(L.八x(3);編寫程序(test2.m)L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160

11、 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;data=K;L;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224256 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata

12、)解得答案為：x =1.2246S4612a = 1.2246, = 0.4612,代.1276即 匕-0.1276根據(jù)程序：L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;K,L=mes

13、hgrid(K,L);Y=1.2246.*(K.A(0.4612).*(L.A(-0.1276);mesh(K,L,Y);做出圖像為：六、模型評(píng)價(jià)根據(jù)程序：L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458

14、458 454;t=1900:1926;data=K;L;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,+,t,qdata1,o)做出擬合值對(duì)比圖像如下：2.6-+2.42.2+ + + +1.81.6+ +41.41.219001905

15、19101915192019251930七、模型改進(jìn)即推廣經(jīng)過(guò)以上計(jì)算出的 為小于0的數(shù)，但是從對(duì)比圖像可以看出擬合效果并不理想。這是因?yàn)椴煌瑲v史時(shí)期 科技發(fā)展?fàn)顟B(tài)不同，而不同科技發(fā)展水平條件下的道格拉斯函數(shù)的系數(shù)是不同的，因此我們擬合的函數(shù)需 要根據(jù)不同時(shí)期科技發(fā)展水平整合。而在不知曉科技發(fā)展水平的情況下，只能通過(guò)已知數(shù)據(jù)分析，因此我 做出下圖對(duì)比Q、K、L增長(zhǎng)情況。t=1900:1926;L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182

16、 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;plot(t,K)hold onplot(t,L)hold onplot(t,qdata)知大致可以把數(shù)據(jù)

17、分為3段，即1900-1906,然后是1907-1910,1911-1926,3段，再次做擬合:編寫M文件curveu.m function q=curvev(x,data) K=data(1,:);L=data(2,:);q=x(1).*(K.八x(2).*(L.八x(3);編寫程序(test2.m)i . L=1e-2*105 108 118 122 117 130 139;K=1e-2*104 106 116 122 127 137 144; data=K;L;qdata=1e-2*105 118 129 130 130 142 150; x0=0.2,0.5,0.5;x=lsqcurv

18、efit(curvev,x0,data,qdata) 解之得：x 1=1.06100.31010.7404ii.L=1e-2* 147 131 143 158;K=1e-2* 153 157 205 251;data=K;L;qdata=1e-2* 152 146 160 169;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata) 解之得：x 2=1.24240.17950.3215iii.L=1e-2* 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2

19、* 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454; data=K;L;qdata=1e-2* 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258; x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之的：x3 =1.48000.4116-0.3425又由上x(chóng)1,x2，已滿足條件，而x3不滿足，則可對(duì)x3段再進(jìn)行分段。則觀察1911-1926圖像如下：t=1911:1926;L=1e-2*

20、 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2* 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454; qdata=1e-2* 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;plot(t,K)hold onplot(t,L)hold onplot(t,qdata)知需將 1910-1926 再分為 1911-1920,1921-1923,192

21、4-1926三段，再次擬合:L=1e-2* 159 166 168 165 162 186 193 196 195 190;K=1e-2* 263 274 282 324 324 361 410 436 477 475;data=K;L;qdata=1e-2* 181 193 195 201 200 209 196 220 212 216;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之得：x 31=1.53370.21700.0023又因?yàn)槎啻螌?shí)驗(yàn)，1920-1926這段時(shí)間內(nèi)的點(diǎn)并不一定滿足所求的道格拉斯公式，因此將其視為奇異點(diǎn)。L=

22、1e-2* 158 167 182;K=1e-2* 454 454 458;data=K;L;qdata=1e-2* 208 224 256;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)L=1e-2* 160 161 164;K=1e-2* 458 458 454;data=K;L;qdata=1e-2* 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之得：x32 =0.01383.18210.6371則綜上可分為3段函數(shù)擬合，所求出的a,a, P如下:a =

23、 1.0610a1 = 0.31011 (1900 -1906)c1p = 0.740412 =a2P2a3a3p3T.2424=0.1795(1907 -1910)=0.3215=1.5337=0.2170 ”911-1920)=0.0023因?yàn)?919以后因?yàn)榭萍及l(fā)展加快則道格拉斯模型不再適用，因此則會(huì)有1919之后的點(diǎn)為奇異點(diǎn)，與道格 拉斯模型擬合并不符合。則再次根據(jù)下列程序t=1900:1906L=1e-2*105 108 118 122 117 130 139;K=1e-2*104 106 116 122 127 137 144;data=K;L;qdata=1e-2*105 118 129 130 130 142 150;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,+,t,qdata1,o)hold ont=1907:1910L=1e-2* 147 131 143 158;K=1e-2* 153 157 205 251;data=K;L;qdata=1e-2* 152 146 160 169;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,d