实验9-人口预测与数据拟合

1、Mathematics Laboratory阮小娥博士Experiments in Mathematics刘康民数学实验数学实验办公地址：理科楼办公地址：理科楼214实验实验13 13 人口预测与数据拟合人口预测与数据拟合2 2、了解利用最小二乘法进行、了解利用最小二乘法进行数据拟合的基本思想，掌握数据拟合的基本思想，掌握用数据拟合法寻找最佳拟合用数据拟合法寻找最佳拟合曲线的方法。曲线的方法。3、了解多元函数的极值在、了解多元函数的极值在数据拟合法中的应用。数据拟合法中的应用。实验目的实验目的1 1、学会用、学会用MATLABMATLAB软件进行数软件进行数据拟合。据拟合。4、通过对实际问题进

2、行分、通过对实际问题进行分析研究，初步掌握建立数析研究，初步掌握建立数据拟合数学模型的方法。据拟合数学模型的方法。据人口统计年鉴，知我国从据人口统计年鉴，知我国从19491949年至年至19941994年人口数据资料如下：年人口数据资料如下： ( (人口数单位为：百万人口数单位为：百万) )（1 1）在直角坐标系上作出人口数的图象。）在直角坐标系上作出人口数的图象。（2 2）建立人口数与年份的函数关系，并估算）建立人口数与年份的函数关系，并估算19991999年年的人口数。的人口数。实验问题实验问题年份年份19491954 1959 1964 1969人口数人口数 541.67602.66 6

3、72.09 704.99 806.71 年份年份 1974 1979 1984 1989 1994人口数人口数 908.59 975.42 1034.751106.761176.74 如何确定如何确定a,b?baxy 线性模型线性模型1 、曲线拟合、曲线拟合问题的提法问题的提法: 已知一组（二维）数据，即平面上的已知一组（二维）数据，即平面上的n 个点个点),(iiyx， ixni, 2 , 1L互不相同，寻求一个函数（曲线）互不相同，寻求一个函数（曲线）)(xfy ，使使)(xf在某种准则下与所有数据点最为接近，即曲线拟合在某种准则下与所有数据点最为接近，即曲线拟合得最好，如图得最好，如图:

4、 xy0+i ),(iiyx)(xfy 一、曲线拟合一、曲线拟合 niiiniiyxf1212)( 确定确定f(x)使得使得 达到最小达到最小 最小二乘准则最小二乘准则 数据插值数据插值已知一组（二维）数据，即平面上的已知一组（二维）数据，即平面上的n 个点个点),(iiyx， ixni, 2 , 1L互不相同，寻求一个函数（曲线）互不相同，寻求一个函数（曲线）)(xfy , 1,)(,)(1niyxfxxxfiinLL 处满足处满足在观察点在观察点使使数据插值数据插值、 用什么样的曲线拟合已知数据用什么样的曲线拟合已知数据?拟合基函数：拟合基函数：画图观察画图观察； 理论分析理论分析12,m

5、xxxeee 指数函数类：指数函数类： 1, ,mxx幂函数类：幂函数类： 12( ),( ),( )mxxx拟合函数：拟合函数：1122( )( )( )( )mmf xaxaxaxLL多项式拟合多项式拟合指数函数指数函数拟合拟合三角函数类：三角函数类： 1,sin ,cos ,sin,cosxxmxmx拟合函数确定后，问题转化为最小值问题：拟合函数确定后，问题转化为最小值问题： 2111min()()nkmmkkkaxaxy LL、 拟合函数组中系数的确定拟合函数组中系数的确定达到最小。达到最小。 niiiniiyaxaaaJ12211221),( 为为此此，只只需需利利用用极极值值的的必

6、必要要条条件件)2 , 1(0 kaJk得得到到 关关于于21,aa的的线线性性方方程程组组， 21,aa,)(21为例为例以以axaxf 12aa即即确确定定 ，使使得得 0) 20) 2121121niiiiniiiyaxaxyaxaMATLAB中的数据拟合命令简单使用格式：简单使用格式： X=lsqcurvefir(fun,x0,xdata,ydata); X,res=lsqcurvefir(fun,x0,xdata,ydata);功能：对数据功能：对数据xdata,ydata按函数按函数fun(文件文件)，以以x0为初值做最小二乘拟合。返回为初值做最小二乘拟合。返回fun中的系中的系数

7、向量和残差平方和的范数。数向量和残差平方和的范数。曲线拟合：曲线拟合：lsqcurvefitlsqcurvefit曲线拟合举例曲线拟合举例例：已知观测数据如下表所示求最小二乘意义的曲线f(x)拟合，其中：X00.10.20.30.40.50.60.70.80.91.0Y3.13.273.814.55.1867.058.569.6911.2513.1723( )xf xaebxcx实验过程：u编写拟合函数文件：编写拟合函数文件：function f=nihefun(x,xdata)function f=nihefun(x,xdata)f=x(1)f=x(1)* *exp(xdata)+x(2)e

8、xp(xdata)+x(2)* *xdata.2+x(3)xdata.2+x(3)* *xdata.3xdata.3u保存为保存为nihefun.mnihefun.mu编写下面的程序：编写下面的程序：sy09_01.msy09_01.m xdata=0:0.1:1; xdata=0:0.1:1; ydata= ydata=3.1 3.27 3.81 4.5 5.18 6 7.05 8.56 9.69 11.25 13.17; x0=0 0 0; x,resn=lsqcurvefit(nihefun,x0,xdata,ydata) xdata=0:0.01:1;nihef=x(1)*exp(xd

9、ata)+x(2)*xdata.2+x(3)*xdata.3; pause plot(xdata,nihef) 00.20.40.60.812468101214MATLAB中的数据拟合命令简单使用格式：简单使用格式： p=polyfit(xdata,ydata,m); 功能：对功能：对xdata,ydata用用m次多项式拟合。次多项式拟合。 返回：拟合多项式的系数向量返回：拟合多项式的系数向量p。（从高到低）。（从高到低） y0=polyval(p,x0)可以求得多项式在可以求得多项式在x0处的处的值。值。 由于高次多项式曲线变化不稳定，所以多项由于高次多项式曲线变化不稳定，所以多项式次数的选

10、取不宜过高。式次数的选取不宜过高。 多项式拟合：多项式拟合：polyfitpolyfit多项式拟合举例多项式拟合举例例：已知观测数据如下表所示分别用3次和6次多项式拟合这些数据点。X00.10.20.30.40.50.60.70.80.91.0Y-.4471.9783.286.167.087.347.669.569.489.311.2实验过程：编写下面的程序：编写下面的程序：sy09_02.msy09_02.m x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2; plot(x,y,k.,markersiz

11、e,25); p3=polyfit(x,y,3); p6=polyfit(x,y,6); t=0:0.01:1.2; s=polyval(p3,t); s1=polyval(p6,t); hold on plot(t,s,k-,linewidth,2) plot(t,s1,k-,linewidth,2) grid on 二、人口预测线性模型二、人口预测线性模型对于开始提出的实验问题，代入数据，计算得对于开始提出的实验问题，代入数据，计算得1215,27754aa 从而得到人口数与年份的函数关系为从而得到人口数与年份的函数关系为2775415 xy把把x=1999代如，估算出代如，估算出1999

12、年的人口数为年的人口数为 y=1252.1（百万）（百万）12.52亿亿1999年实际人口数量为年实际人口数量为.亿。亿。线性预测模型线性预测模型 英国人口学家英国人口学家MalthusMalthus根据百余年的人口统计资根据百余年的人口统计资料，于料，于17981798年提出了著名的人口自然增长的年提出了著名的人口自然增长的指数增指数增长模型长模型。 btaetx 三、人口预测的三、人口预测的Malthus模型模型 0)0(,xxrxdtdxtrextx0)( 基本假设基本假设 : : 人口人口( (相对相对) )增长率增长率 r 是常数是常数x(t) 时刻时刻t 的的人口人口,t=0时人口

13、数为时人口数为x0 btatx ln指数增长模型指数增长模型实际中，常用实际中，常用1. 1. 由前由前100100年的数据求出美国的人口增长年的数据求出美国的人口增长Malthus模型模型。2. 2. 预测后预测后100100年（每隔年（每隔1010年）的人口状况。年）的人口状况。3. 3. 根据预测的人口状况和实际的人口数量根据预测的人口状况和实际的人口数量, ,讨论人讨论人口模型的改进情况。口模型的改进情况。美国美国17901790年年19801980年每隔年每隔1010年的人口记录年的人口记录226.5204.0179.3150.7131.7123.2106.592.076.062.9

14、人口人口( (百万百万) )1980197019601950194019301920191019001890年份年份50.238.631.423.217.112.99.67.25.33.9人口人口( (百万百万) )1880187018601850184018301820181018001790年份年份例例解：解： btaetx 设设btax )ln(2101)ln(),(iiixbtabaJ 取得最小值取得最小值. .其中其中, , iitxx 表示人口数量表示人口数量。it表示年份表示年份, ,解方程组解方程组: 0ln210ln10211011012101101101iiiiiiiiii

15、itxbtatbJxbtaaJ即得参数即得参数ba,的值的值.ba,使得使得问题转化为求参数问题转化为求参数 sy09_3.msy09_3.m % % % This program is to predict the number of This program is to predict the number of population %population %format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;187t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;

16、1880;0;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;197t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;0;1980;x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;2;x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.x2=62.9;76.0;92.0;106.5;123.2;131.7;15

17、0.7;179.3;204.0;226.5;3;204.0;226.5;lnx1=log(x1); lnx2=log(x2);lnx1=log(x1); lnx2=log(x2);a12=sum(t1);a11=10;a21=a12;a22=sum(t1.2); d1=sum(lnx1);d2=sum(lnx1.*t1);A=a11,a12;a21,a22;D=d1;d2;ab=inv(A)*D;disp(a=);disp(ab(1);disp(b=);disp(ab(2);forfor i=1:10 i=1:10 xx1(i)=exp(ab(1)+ab(2) xx1(i)=exp(ab(1)

18、+ab(2)* *t1(i);t1(i);endendforfor i=1:10 i=1:10 xx2(i)=exp(ab(1)+ab(2) xx2(i)=exp(ab(1)+ab(2)* *t2(i);t2(i);endendplot(t1,x1,rplot(t1,x1,r* *-,t1,xx1,b+-, -,t1,xx1,b+-, t2,x2,gt2,x2,g* *-,t2,xx2,m+-);-,t2,xx2,m+-);a= -49.79535457790735b=0.02859807120038拟合结果拟合结果表明：表明： 人口的人口的指数模指数模型在短期内基型在短期内基本上能比较准本上

19、能比较准确地反映人口确地反映人口自然增长的规自然增长的规律，但长期预律，但长期预测误差很大，测误差很大，需要修正预测需要修正预测模型。模型。拟合曲线拟合曲线原始数据曲线原始数据曲线四、人口预测的四、人口预测的Logistic模型模型如果人口的增长符合如果人口的增长符合Malthus模型，模型，则则人口增长到一定数量后，增长率下降的原因：人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用资源、环境等因素对人口增长的阻滞作用 1838 1838年，荷兰生物学家年，荷兰生物学家VerhulstVerhulst对对MalthusMalthus模型作模型作了进一步分析后指出：了

20、进一步分析后指出：导致上述不符合实际情况的导致上述不符合实际情况的主要原因是未能考虑主要原因是未能考虑“密度制约密度制约”因素。因素。( )tx t 即最终导致地球上人口爆炸，这与实际是不相符的。即最终导致地球上人口爆炸，这与实际是不相符的。且阻滞作用随人口数量增加而变大且阻滞作用随人口数量增加而变大r是是x的减函数的减函数假设假设r固有增长率固有增长率(x很小时很小时)k人口容量（资源、环境能容纳的最大数量）人口容量（资源、环境能容纳的最大数量）)1()(kxrxr krs 0)( kr)1()(mxxrxxxrdtdx 000( )()rtkxx tkx ex 例2 中国人口预测year

21、popu popu_fit 1971 8.523 8.619 1972 8.718 8.7601973 8.921 8.9011974 9.086 9.0421975 9.242 9.183 1976 9.372 9.324 1977 9.497 9.466 1978 9.626 9.607 1979 9.754 9.749 1980 9.871 9.8901981 10.007 10.031 1982 10.165 10.172 year popu popu_fit 1983 10.301 10.3131984 10.436 10.4531985 10.585 10.5931986 10.7

22、51 10.7321987 10.930 10.8711988 11.103 11.0091989 11.270 11.1471990 11.433 11.284中国人口Logistic模型zgrk_lgstic.mclear;clc;clfload zgrk_data;k=19;cpopu1=1./cpopu-1/k;n=15;a=sum(year(1:n);b=sum(year(1:n).*year(1:n);c=sum(log(cpopu1(1:n);d=sum(year(1:n).*log(cpopu1(1:n);A=n a;a b;B=c;d;p=inv(A)*Bx=1971:200

23、6;y=1./(1/k+exp(p(1)+p(2)*x);cpopu_new=11.582,11.717,11.852,11.985,12.112,12.239,. 12.363,12.476,12.519,12.634,12.763,12.845,12.923,. 12.999,13.076,13.145;cpoputj=cpopu,cpopu_new; scatter(x,cpoputj,k*)axis(1971 2010 6 15);grid;hold on；plot(x,y,r-,linewidth,2)fprintf( year popu popu_fitn);for i=1:20

24、year=i+1970; fprintf(%6.0f %8.3f %8.3fn,year,cpoputj(i),y(i)endfprintf(=n);fprintf( year opu_tj popu_predn);for i=21:36 year=i+1970; fprintf(%6.0f %8.3f %8.3fn,year,cpoputj(i),y(i);end中国人口预测的中国人口预测的Logistic模型模型结果：假设中国人口容纳量为19亿时，Logistic模型为：155.9493 0.02981( )19etN t 年年 份份201020152020人口数（亿）13.798914.

25、342214.8458多项式拟合人口模型多项式拟合人口模型sy09_04.m% % This program is to predict the model of population by This program is to predict the model of population by 4-degree polynomial% 4-degree polynomial% % %prog42.m%prog42.m%format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t1=1790;1800

26、;1810;1820;1830;1840;1850;1860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t=t1;t2;t=t1;t2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P2=62.9;76.0;92.0;106.5;123.2;131.7;150

27、.7;179.3;204.0;2P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;26.5;P=P1;P2;P=P1;P2;n=4; n=4; % The degree of the fitting polynomial% The degree of the fitting polynomial%a,s=polyfit(t1,P1,n);a,s=polyfit(t1,P1,n);y=polyval(a,t);y=polyval(a,t);% % a a is the coefficients vector from n-de

28、gree to 0- is the coefficients vector from n-degree to 0-degree%degree%plot(t,P,rplot(t,P,r* *-,t,y,b+-);-,t,y,b+-);a =1.0e+006 * -0.00000000000014 0.00000000107892 -0.00000304878595 0.00381927346813 -1.79012132225427仿真结果表明仿真结果表明, , 人口增加的模型用多项式拟合能人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预比较准确地反映人口自然增长的规律，对

29、长期预测具有指导意义。测具有指导意义。例例3: 3: 海底光缆线长度预测模型海底光缆线长度预测模型某一通信公司在一次某一通信公司在一次施工中施工中, ,需要在水面宽需要在水面宽为为2020m m的河沟底沿直线的河沟底沿直线走向铺设一条沟底光走向铺设一条沟底光缆缆. .在铺设光缆之前需在铺设光缆之前需要对沟底的地形做初要对沟底的地形做初Bix2468101214161820986420ADCih0 x20 x探测到一组等分点位置的深度数据如下表所示探测到一组等分点位置的深度数据如下表所示. .25步探测步探测, ,从而估计所需光缆的长度从而估计所需光缆的长度, ,为工程预算为工程预算提供依据提供

30、依据. .基本情况如图所示基本情况如图所示. .10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度深度(m)9876543210分点分点2121个等分点处的深度个等分点处的深度(1) (1) 预测通过这条河沟所需光缆长度的近似值预测通过这条河沟所需光缆长度的近似值. .(2) (2) 作出铺设沟底光缆的曲线图作出铺设沟底光缆的曲线图. .024681012141618207891011121314解：解： 用用12次多项式函数拟合光缆走势的曲线图如下次多项式函数拟合光缆走势的曲线图如下仿真结果表仿真结果表明明,拟合曲线拟合曲线能较准确地能较准确