实验十三：数据拟合与数据差值
练习三
1.对示例10中给出的人口统计数据,用不同次数的多项式拟合,预测相应5个年份的人口数,你会发现什么?
clc;clear;
x0=1971:1990;
y0=[8.523,8.718,8.921,9.086,9.242,9.372,9.497,9.626,9.754,9.871,10.007,10.165,10.301,10.436,10.585,10.751,10.930,11.103,11.270,11.433];
p1=polyfit(x0,y0,2);
x=1995:5:2015;
y1=polyval(p1,x)
p2=polyfit(x0,y0,5);
y2=polyval(p2,x)
y1 =
12.1797 12.9893 13.8224 14.6789 15.5587y2 =
y2 =
12.0502 11.6070 8.7101 1.3216 -13.3177
当次数较小时，预测还较准确；当次数稍大时，预测不合理。
2.1790年到1980年间美国人口数的统计数据如表13.13所示
表13.13 美国人口统计数据
| 年份 | 1790 | 1800 | 1810 | 1820 | 1830 | 1840 | 1850 | 1860 | 1870 | 1880 |
| 人口数/百万 | 3.9 | 5.3 | 7.2 | 9.6 | 12.9 | 17.1 | 23.2 | 31.4 | 38.6 | 50.2 |
| 年份 | 1890 | 1900 | 1910 | 1920 | 1930 | 1940 | 1950 | 1960 | 1970 | 1980 |
| 人口数/百万 | 62.9 | 76 | 92.0 | 105.7 | 122.8 | 131.7 | 150.7 | 179.3 | 203.2 | 226.5 |
(1)根据表13.13中的数据,分别用不同次数多项式拟合美国人口数增长的近似曲线;
clc;clear;
x0=1790:10:1980;
y0=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92.0,105.7,122.8,131.7,150.7,179.3,203.2,226.5];
x=1790:1980;
p1=polyfit(x0,y0,2);
y1=polyval(p1,x);
figure(1)
plot(x0,y0,'.',x,y1);
p2=polyfit(x0,y0,10);
y2=polyval(p2,x);
figure(2)
plot(x0,y0,'.',x,y2);
(2)根据表13.13中的数据,建立符合马尔萨斯模型的美国人口数增长模型;
plot(x0,y0,'.',x,y2);
%求a和b
a0=[0,0];
f1=@(a,t)exp(a(1)+a(2)*t);
a=nlinfit(x0,y0,f1,a0);
xx=1790:1980;
g1=exp(a(1)+a(2)*xx);
plot(x0,y0,'.',xx,g1);
(3)设美国人口总体容纳量为4.5亿,试用逻辑斯谛模型建立美国人口增长模型;
k=450;
%求a和b
a0=[0,0];
f2=@(b,t)1./(k^-1+exp(-b(1)-b(2)*t));
b=nlinfit(x0,y0,f2,a0);
xx=1790:1980;
g2=1./(k^-1+exp(-b(1)-b(2)*xx));
plot(x0,y0,'.',xx,g2);
(4)分别用上述三种方法预测2000年,2005年,2010年,2015年,2020年美国的人口数,并对不同方法的预测结果进行比较分析.
ycx=2000:5:2020;
ycy1=polyval(p1,ycx)
ycy2=polyval(p2,ycx)
g11=@(x)exp(a(1)+a(2)*x);
ycy3=g11(ycx)
g22=@(x)1./(k^-1+exp(-b(1)-b(2)*x));
ycy4=g22(ycx)
ycy1 =
1.0e+02 *
2.739705263157921 2.873557792207830 3.010662884483936 3.151020539986312 3.294630758714993
ycy2 =
1.0e+03 *
0.059429687500000 -0.103941406250000 -0.364816406250000 -0.761148437500000 -1.342269531250000
ycy3 =
1.0e+02 *
3.252308701187482 3.513505561833679 3.795679459495618 4.100515085485712 4.429832430721169
ycy4 =
1.0e+02 *
2.741520181295475 2.853817611979104 2.962964306661279 3.068504260808794 3.170050165270613
2次多项式拟合
| 年份 | 2000 | 2005 | 2010 | 2015 | 2020 |
| 人口数/百万 | 273.97 | 287.36 | 301.07 | 315.10 | 329.46 |
10次多项式拟合
| 年份 | 2000 | 2005 | 2010 | 2015 | 2020 |
| 人口数/百万 | 59.43 | -103.94 | -364.82 | -761.15 | -1342.27 |
马尔萨斯模型模型
| 年份 | 2000 | 2005 | 2010 | 2015 | 2020 |
| 人口数/百万 | 325.23 | 351.35 | 379.57 | 410.05 | 442.98 |
逻辑斯蒂模型
| 年份 | 2000 | 2005 | 2010 | 2015 | 2020 |
| 人口数/百万 | 274.15 | 285.38 | 296.30 | 306.85 | 317.01 |
3.由于地球围绕太阳公转,所以气温具有周期性,某地气象台测得当地月平均气温如表13.14所示
表13.14 某地全年月平均气温数据
| 月份 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
| 平均气温/℃ | 3.1 | 3.8 | 6.9 | 12.7 | 16.8 | 20.5 | 24.5 | 25.9 | 22 | 16.1 | 10.7 | 5.4 |
试建立该地全年温度与时间的经验公式.(提示:可考虑用傅里叶级数.)
这里我们说明一下Curve Fitting Tool 的使用方法（简称为cftool）
首先我们应该如何安装cftool工具箱呢？
在matlab选项栏中点击APPà获取更多APP，之后出来一个名叫附加功能资源管理器的窗口，在搜索框中输入Curve Fitting Tool，选择该工具箱安装即可。（如果用户没有登录的话需要登录账户）
安装完成之后在命令行窗口输入cftool，出现一个曲线拟合器的窗口（如下图）
点击上方的选择数据选项。（点击之前需要在原有文件中运行一下程序）
clc;clear;
month=1:12;
t=[3.1,3.8,6.9,12.7,16.8,20.5,24.5,25.9,22,16.1,10.7,5.4];
在选择拟合数据窗口中填写相应的信息即可关闭。
针对本题而言，我们选择傅里叶级数拟合该离散点。
首先，我们完成以上操作，点击上方傅里叶选项即可出现运行结果。
除此之外，我们还可以生成相应的代码，导出图像或者导出到工作区等等。
导出代码：
function [fitresult, gof] = createFit(month, t)
%CREATEFIT(MONTH,T)
% 创建一个拟合。
%
% 要进行 '傅里叶拟合' 拟合的数据:
% X 输入: month
% Y 输出: t
% 输出:
% fitresult: 表示拟合的拟合对象。
% gof: 带有拟合优度信息的结构体。
%
% 另请参阅 FIT, CFIT, SFIT.
% 由 MATLAB 于 16-Feb-2024 19:23:55 自动生成
%% 拟合: '傅里叶拟合'。
[xData, yData] = prepareCurveData( month, t );
% 设置 fittype 和选项。
ft = fittype( 'fourier1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0 0 0 0.571198664289053];
% 对数据进行模型拟合。
[fitresult, gof] = fit( xData, yData, ft, opts );
% 绘制数据拟合图。
figure( 'Name', '傅里叶拟合' );
h = plot( fitresult, xData, yData );
legend( h, 't vs. month', '傅里叶拟合', 'Location', 'NorthEast', 'Interpreter', 'none' );
% 为坐标区加标签
xlabel( 'month', 'Interpreter', 'none' );
ylabel( 't', 'Interpreter', 'none' );
grid on
下图是相应的数据结果：
4.1601年,德国天文学家开普勒(kepler)发表了行星运行第三定律: T=Cx3/2
,其中T为行星绕太阳旋转一周的时间(单位:天),x表示行星到太阳的平均距离(单位:106
km),并测得水星、金星、地球、火星的(x,T)数据分别为(58,88),(108,225),(150,365),(228,687).
(1)用最小二乘法估计C的值;
clc;clear;
x=[58,108,150,228];
y=[88,225,365,687];
f=inline('a*x*(3/2)','a','x');
a0=0;
[a,r]=lsqcurvefit(f,a0,x,y)
a =
1.7934
r =
1.5831e+04
(2)分别作出上述数据点的直线、抛物线、3次多项式、4次多项式拟合,求出残差平方和,并比较优劣;
注：此题应该是题目有问题：无法进行四次拟合，因为题目中只给了四个数据。
我们可以先用正常代码试一下
clc;clear;
x=[58,108,150,228];
y=[88,225,365,687];
n=1;%多项式拟合次数，可以修改
p1=polyfit(x,y,n);
xx=58:228;
y1=polyval(p1,xx);
y11=polyval(p1,x);
plot(x,y,'.',xx,y1);
r=0;
for i=1:length(x)
r=r+(y11(i)-y(i))^2;
end
r
上面的代码我只试用了直线拟合。
我们可以试一下刚才用的那个曲线拟合器。
拟合数据：
拟合名称: 一次拟合
多项式 曲线拟合(poly1)
f(x) = p1*x + p2
系数和 95% 置信边界
值 下限 上限
p1 3.5516 2.4921 4.6112
p2 -141.7742 -300.2771 16.7287
拟合优度
值
SSE 1.8833e+03
R 方 0.9905
DFE 2
调整 R 方 0.9857
RMSE 30.6865
其中SSE是残差平方和。
其他的我就不一一复制了。
(3)用函数y=aex+bx+c
对数据点进行曲线拟合,并求出残差平方和.
此题我感觉出的不是很对，因为这个拟合十分奇怪，对初值要求较高；（exp（x）感觉太猛了，个人认为把他改成exp（-x）较好）
clc;clear;
x=[58,108,150,228];
y=[88,225,365,687];
f=inline('a(1)*exp(x)+a(2)*x+a(3)','a','x');
a0=[0,3.5,-141];
[a,r]=lsqcurvefit(f,a0,x,y)
xx=58:228;
g=a(1)*exp(xx)+a(2)*xx+a(3);
plot(x,y,'.',xx,g);
a =
0.0000 3.5000 -141.0000
r =
1181
上边的是按题意来的，故意把初值选成一次多项式拟合的值；
如果我们改成exp（-x）试一下
clc;clear;
x=[58,108,150,228];
y=[88,225,365,687];
f=inline('a(1)*exp(-x)+a(2)*x+a(3)','a','x');
a0=[0,0,0];
[a,r]=lsqcurvefit(f,a0,x,y)
xx=58:228;
g=a(1)*exp(-xx)+a(2)*xx+a(3);
plot(x,y,'.',xx,g);
a =
0 3.5516 -141.7742
r =
1.8833e+03
可以对比得知，此时的a值还有残差平方和和刚才的拟合几乎一样。
推荐下一篇文章：
https://blog.csdn.net/2301_80199493/article/details/136158596?spm=1001.2014.3001.5501本文由作者自创，由于时间原因，难免出现些许错误，还请大家多多指正。创作不易，请大家多多支持。
717
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
