马尔萨斯数据拟合实验

第第#页共12页马尔萨斯数据拟合实验实验内容实验一：数据拟合Malthus人口指数增长模型中的参数1790-1980年间美国每隔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人口（X106）123.2131.7150.7179.3204.0226.5用以上数据检验马尔萨斯（Malthus）人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。试验二：经济增长模型增加生产、发展经济所依靠的主要因素有增加投资、增加劳动力以及技术革新等，在研究国名经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易定量化，作为初步的模型，可认为技术水平不变，只讨论不像资金、劳动力之间的关系。在科学技术发展不快时，如资本主义经济发展的前期,这种模型是有意义的。用Q,K,L分别表示产值、资金、劳动力，要寻求的数量关系Q（K,L）。经过简化假设与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数(*)Q（K丄）=aKaLP,0<°,卩V1(*)式中a,P,a要由紧急统计数据确定。现有美国马萨诸赛州1900-1926年上述三个经济指数的统计数据。如表5.5，试用数据拟合的方法，求出（*）式中的参数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.164.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三、实验环境Windows操作系统；MATLAB7.0.四、实验过程一、问题：实验一二、问题分析因为假设人口增长率为常数，且已知初始时刻的人口数，可根据t时刻人口变化用差商代替微商。当AtT0时，可列出关于t时刻人口变化的微分方程。dx—rx<dt。根据该微分方程求出解，再用数据拟合模型球参数。x(0)=xI0三、假设假设人口增长率为常数，记为r.记时刻t的人口为x(t)(即x(t)为模型的状态变量)。假设初始时刻人口为x0。四、模型建立求解该微分方程得到x(t)的表达式为x£)=x*er*t则可转化为求关于x(t)在t时刻的测量数据来辨识0参数对应的数学模型球函数的最小二乘解。其中a=0,b=x0,则需求r的最小值点。运用最小二函数拟合求出所需的r值。五、模型求解编写M文件curvefun.mfunctionf=curvefun(x,tdata)f=1e+05*39*exp(x*tdata);编写(testl.m)tdata=linspace(0,200,20);cdata=le+05.*[395372961291712323143865026297609201065123213171507179320402265];x0=[0.2];y=lsqcurvefit('curvefun',x0,tdata,cdata)即求出结果为：y=0.0211即r=0.0211根据下述程序作出模型图像：tdata=linspace(0,200,20);cdata=1e+05.*[395372961291712323143865026297609201065123213171507179320402265];plot(tdata,cdata)holdonf=1e+05*39*exp(0.0211*tdata);plot(tdata,f,'r'),title(‘马尔萨斯模型')模型图像：六、模型评价可以看出使用一般曲线拟合，出来走势基本一致，但是还有一定误差，并且由图可以看出该曲线接近一条开口向上的二次曲线。因此可以尝试使用二次多项式对马尔萨斯模型进行拟合。将会大大简化求解的过程。七、模型改进即推广则二次曲线拟合求解如下：编写程序：tdata=linspace(0,200,20);cdata=le+05.*[395372961291712323143865026297609201065123213171507179320402265];p=polyfit(tdata,cdata,2)得到结果：p=1.0e+006*0.0059-0.08265.4634编程作图：tdata=linspace(0,200,20);cdata=1e+05.*[395372961291712323143865026297609201065123213171507179320402265];plot(tdata,cdata)holdonf=1.0e+006*5.4634.*tdata.A2-1.0e+006*0.0826.*tdata+1.0e+006*0.0059；plot(tdata,f,'r'),title(‘马尔萨斯模型1')则可根据上述程序作出二次曲线与点线图如下：可见该二次曲线与马尔萨斯模型基本上是完全重合的。因此儿二次曲线拟合是优于一般的曲线拟合的，对于马尔萨斯模型来说。

一、问题：实验二二、问题分析因为已知生产函数为Q(KL=aK吋,0<?卩<1，且已知初始时刻的相关Q、K、L，可根据上式进行数据拟合，但是需将其转化为2维数组计算即将L、K合成为一个2维数组data再从data中将K、L取出，则可利用一维拟合求出所要的解。三、假设Q,K,L分别表示产值、资金、劳动力认为技术水平不变，只讨论不像资金、劳动力之间的关系四、模型建立五、模型求解编程如下：编写M文件curveu.mfunctionq=curvev(x,data)K=data(1,:);L=data(2,:);q=x(1).*(K.Ax(2)).*(L.Ax(3));编写程序(test2.m)L=le-2*[105108118122117130139147131143158159166168165162186193196195190158167182160161164];K=1e-2*[104106116122127137144153157205251263274282324324361410436477475454454458458458454];data=[K;L];qdata=1e-2*[105118129130130142150152146160169181193195201200209196220212216208224256234245258];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解得答案为：x=1.22460.46121.22460.4612a二1.2246,二0.4612,代.1276即匕-0.1276根据程序：L=1e-2*[105108118122117130139147131143158159166168165162186193196195190158167182160161164];K=1e-2*[104106116122127137144153157205251263274282324324361410436477475454454458458458454];[K,L]=meshgrid(K,L);Y=1.2246.*(K.a(0.4612)).*(L.a(-0.1276));mesh(K,L,Y);做出图像为：

45324532六、模型评价根据程序：L=le-2*[105108118122117130139147131143158159166168165162186193196195190158167182160161164];K=1e-2*[104106116122127137144153157205251263274282324324361410436477475454454458458458454];t=1900:1926;data=[K;L];qdata=1e-2*[105118129130130142150152146160169181193195201200209196220212216208224256234245258];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,'+',t,qdata1,'o')做出拟合值对比图像如下：

2.6-+2.4-+2.2-+1.81.6-+1.41.2「900++419051910■+-++亠19152.6-+2.4-+2.2-+1.81.6-+1.41.2「900++419051910■+-++亠1915192019251930七、模型改进即推广经过以上计算出的为小于0的数，但是从对比图像可以看出拟合效果并不理想。这是因为不同历史时期科技发展状态不同，而不同科技发展水平条件下的道格拉斯函数的系数是不同的，因此我们拟合的函数需要根据不同时期科技发展水平整合。而在不知晓科技发展水平的情况下，只能通过已知数据分析，因此我做出下图对比Q、K、L增长情况。t=1900:1926;L=le-2*[105108118122117130139147131143158159166168165162186193196195190158167182160161164];K=1e-2*[104106116122127137144153157205251263274282324324361410436477475454454458458458454];qdata=1e-2*[105118129130130142150152146160169181193195201200209196220212216208224256234245258];plot(t,K)holdonplot(t,L)holdonplot(t,qdata)知大致可以把数据分为3段，即1900-1906，然后是1907-1910,1911-1926,3段，再次做拟合:编写M文件curveu.mfunctionq=curvev(x,data)K=data(1,:);L=data(2,:);q=x(1).*(K.Ax(2)).*(L.Ax(3));编写程序(test2.m)i.L=1e-2*[105108118122117130139];K=1e-2*[104106116122127137144];data=[K;L];qdata=1e-2*[105118129130130142150];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解之得：x1=1.06100.31010.7404ii.L=1e-2*[147131143158];K=1e-2*[153157205251];data=[K;L];qdata=1e-2*[152146160169];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解之得：x2=1.24240.17950.3215iii.L=le-2*[159166168165162186193196195190158167182160161164];K=le-2*[263274282324324361410436477475454454458458458454];data=[K;L];qdata=1e-2*[181193195201200209196220212216208224256234245258];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解之的：x3=1.48000.4116-0.3425又由上x1,x2，已满足条件，而x3不满足，则可对x3段再进行分段。则观察1911-1926图像如下：t=1911:1926;L=1e-2*[159166168165162186193196195190158167182160161164];K=1e-2*[263274282324324361410436477475454454458458458454];qdata=1e-2*[181193195201200209196220212216208224256234245258];plot(t,K)holdonplot(t,L)holdonplot(t,qdata)知需将1910-1926再分为1911-1920,1921-1923,1924-1926三段，再次拟合:L=1e-2*[159166168165162186193196195190];K=le-2*[263274282324324361410436477475];data=[K;L];qdata=1e-2*[181193195201200209196220212216];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解之得：x31=1.53370.21700.0023又因为多次实验，1920-1926这段时间内的点并不一定满足所求的道格拉斯公式，因此将其视为奇异点。L=1e-2*[158167182];K=1e-2*[454454458];data=[K;L];qdata=1e-2*[208224256];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)L=1e-2*[160161164];K=1e-2*[458458454];data=[K;L];qdata=1e-2*[234245258];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)解之得：x32=0.01383.18210.6371则综上可分为3段函数拟合，所求出的a,a,P如下:a=1.0610a=0.3101)(1900-1906)1P=0.74041a=1.2424=<a=0.1795[(1907-1910)2P=0.32152a=1.5337=0.2170)(1911-1920)3P=0.00233丿因为1919以后因为科技发展加快则道格拉斯模型不再适用，因此则会有1919之后的点为奇异点，与道格拉斯模型拟合并不符合。则再次根据下列程序t=1900:1906L=1e-2*[105108118122117130139];K=1e-2*[104106116122127137144];data=[K;L];qdata=le-2*[105118129130130142150];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,'+',t,qdata1,'o')holdont=1907:1910L=1e-2*[147131143158];K=1e-2*[153157205251];data=[K;L];qdata=1e-2*[152146160169];x0=[0.2,0.5,0.5];x=lsqcurvefit('curvev',x0,data,qdata)qdata1=curveu(x,dat