数据拟合

%%%%%%%
数据拟合
根据一组二维数据，即平面上的若干点，要求确定一个一元函数y =
f(x)，即曲线，使这些点与曲线总体来说尽量接近。

这就是数据拟合成曲线的思想，简称为曲线拟合(fitting a curve)。

曲线拟合其目的是根据实验获得的数据去建立因变量与自变量之间有效的经验函数关系，为进一步的深入研究提供线索。

本章的目的，掌握一些曲线拟合的基本方法，弄清楚曲线拟合与插值方法之间的区别,学会使用MATLAB软件进行曲线拟合。

§5.1 引例
拟合问题引例一电阻问题
已知热敏电阻电阻值与温度的数据：
求温度为63度时的电阻值。

拟合问题引例二给药问题
一种新药用于临床之前，必须设计给药方案。

药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。

一室模型：将整个机体看作一个房室，称中心室，室内血药浓度是均匀的。

快速静脉注射后，浓度立即上升；然后迅速下降。

当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。

临床上，每种药物有一个最小有效浓度c 1和一个最大有效浓度c 2。

设计给药方案时，要使血药浓度 保持在c 1~c 2之间。

本题设c 1=10，c 2=25(ug/ml).
要设计给药方案,必须知道给药后血药浓度随时间变化的规律。

从实验和理论两方面着手：
在实验方面, t=0时对某人用快速静脉注射方式一次注入该药物300mg 后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表:
1. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。

2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。

§5.2 最小二乘法
给定平面上的点(x i , y i )，（i = 1,2,…,n ），进行曲线拟合有多种方法，其中最小二乘法是解决曲线拟合最常用的方法。

最小二乘法的原理是：求 ∑∑==-==n
i i i n
i i y x f x f 1
21
2])([),(δδ使 达到最小。

如图1所示，其中δi 为点(x i ,y i )与曲线y=f (x )的距离。

曲线拟合的实际含义是寻求一个函数y=f (x )，使f (x )在某种准则下与所有数据点最为接近，即曲线拟合得最好。

最小二乘准则就是使所有散点到曲线的距离平方和最小。

拟合时选用一定的拟合函数f (x ) 形式，设拟合函数可由一些简单的“基函数”（例如幂函数，三角函数等等） )(),...,(),(10x x x m ϕϕϕ来线性表示：
)(...)()()(1100x c x c x c x f m m ϕϕϕ+++=
图1 曲线拟合示意图
现在要确定系数c 0，c 1，…，c m ，使δ达到极小。

为此，将f (x )的表达式代入δ中，δ就成为c 0，c 1，…，c m 的函数，求δ的极小，就可令δ对 c i 的偏导数等于零，于是得到m +1个方程组，从中求解出c i 。

通常取基函数为1，x ，x 2，x 3，…，x m ，这时拟合函数f (x )为多项式函数。

当m =1时，f (x ) = a + bx ，称为一元线性拟合函数，它是曲线拟合最简单的形式。

除此之外，常用的一元曲线拟合函数还有双曲线f (x ) = a + b/x ，指数曲线f (x ) = a e b x 等，对于这些曲线，拟合前须作变量代换，转化为线性函数。

已知一组数据，用什么样的曲线拟合最好呢？可以根据散点图进行直观判断，在此基础上，选择几种曲线分别作拟合，然后比较，观察哪条曲线的最小二乘指标δ最小。

§5.3 曲线拟合的MATLAB 实现
MATLAB 软件提供了基本的曲线拟合函数的命令：
22
f=a +a x
1．多项式函数拟合： a = polyfit （xdata ，ydata ，n ）
其中n 表示多项式的最高阶数，xdata ，ydata 为要拟合的数据，它是用数组的方式输入。

输出参数a 为拟合多项式y = a 1x n + … + a n x + a n +1的系数a = [a 1, …, a n , a n +1]。

多项式在x 处的值y 可用下面程序计算。

y = polyval (a, x)
2．一般的曲线拟合： p = curvefit （‘Fun ’p0，xdata ，
ydata ）
其中Fun 表示函数Fun (p, xdata)的M-文件，p0表示函数的初值。

curvefit 命令的求解问题形式是：
min {p } sum {(Fun (p, xdata)-ydata).^2}
若要求解点x 处的函数值可用程序f = Fun(p, x) 计算。

例如已知函数形式 dx bx ce ae y --+=，并且已知数据点(x i , y i ), i = 1,2,…, n ，要确定四个未知参数a, b, c, d 。

使用curvefit 命令，数据输入xdata = [x 1,x 2, …, x n ]; ydata = [y 1,y 2, …, y n ];初值输入p0 = [a 0,b 0,c 0,d 0]; 并且建立函数
dx bx ce ae y --+=的M-文件（Fun.m ）。

若定义p 1 = a , p 2 = b , p 3 = c , p 4 = d , 则输出p = [p 1, p 2, p 3, p 4]。

3．引例1求解
根据热敏电阻电阻值与温度的数据：
首先作出散点图
t=[20.5 32.7 51.0 73.0 95.7]; r=[765 826 873 942 1032]; plot(t,r)
从散点图看出电阻值和温度之间的关系近似于线性。

因此设R=a1 t +a2 a=polyfit(t,r,1)
a =
3.3987 702.0968
表示R=3.3987 t +702.0968
63度时的电阻值为
r=polyval(a,63)
r =
916.2174
4．引例2求解：
（1）模型假设
1. 机体看作一个房室，室内血药浓度均匀——一室模型
2. 药物排除速率与血药浓度c成正比，比例系数 k(>0)
3. 血液容积 v, t=0时注射剂量 d, 血药浓度即为 d/v.
（2）模型建立 由假设2得：
kt e v
d t c kc dt
dc
-=
⇒-=)( 由假设3得：v d c /)0(=
在此，d=300mg ，t 及c （t ）在某些点处的值见前表，需经拟合求出参数k 、v
}2
/,)
/ln(,,ln )/ln(ln )(12
121a kt
e d v a k a t a y v d a k a c y kt v d c e v
d t c =-=+=⇒
=-==-=⇒=
-
d=300;
t=[0.25 0.5 1 1.5 2 3 4 6 8];
c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; y=log(c);
a=polyfit(t,y,1) k=-a(1)
v=d/exp(a(2)) 计算结果： a =
-0.2347 2.9943 k =
0.2347 v =
15.0219
画出血药浓度随时间的变化规律图 t1=[0:0.1:8];
ct=(d/v)*exp(-k*t1); plot(t,c,'o',t1,ct,'g-')
• 设每次注射剂量D, 间隔时间τ • 血药浓度c(t) 应c 1≤ c(t) ≤ c 2 • 初次剂量D 0 应加大 给药方案记为：{}τ，，D D 0 1.)(,1220c c v D vc D -==
2.
12
21ln 1c c
k e c c k =⇒=-ττ 其中c1=10，c2=25，k=0.2347，v=15.02 计算结果：
9.3,3.225,5.3750===τD D
给药方案：
)(4),(3.225),(5.3750h mg D mg D ===τ
即:
首次注射 375 mg ，
其余每次注射 225 mg，
注射的间隔时间为 4 小时。

5.拟合与插值的比较
问题：给定一批离散的数据点，需确定满足特定要求的曲线或曲面,
从而获取整体的规律。

即通过"窥几斑"来达到"知全豹"。

解决方案：
若要求所求曲线（面）通过所给所有数据点，就是插值问题；
若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。

从几何意义上看，拟合是给定了空间中的一些点，找到一个已知形式的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。

附：书上引例：浓度变化规律的计算
x=[1:16];
y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10. 55,10.58,10.6];
plot(x,y,'o')
p=polyfit(x,y,2)
计算结果：
p = -0.0445 1.0711 4.3252 %二次多项式的系数
从而得到某化合物的浓度y与时间x的拟合函数：
2
x
=
.4x
+
y-
3252
.0
0711
0445
.1
对函数的精度如何检测呢？仍然以图形来检测，将散点与拟合曲线画在一个画面上。

xi=linspace(0,16,160);
yi=polyval(p,xi);
plot(x,y,'o',xi,yi)
（这两段程序在此无法运行，需在MATLAB的命令窗口中运行）
参见图2。

由此看出上述曲线拟合是比较吻合的。

图2 浓度y的拟合曲线与实测数据(o)的比较下面我们来看看曲线拟合与曲线插值有什么区别：
MATLAB程序如下：
t=interp1(x,y,xi,'linear');
plot(x,y,'+',xi,t, 'r:')
hold on
plot(xi,yi)
该程序运行后，得到图3所示的结果。

其中，标有"+"的是已知数据点，连接数据点的虚线是线性插值函数曲线，光滑的函数曲线是最佳拟合曲线。

图3 拟合曲线与线性插值曲线的比较
在MATLAB的NAG Foundation Toolbox中也有一些曲面拟合函数，如
e02daf，e02cf，e02def可分别求出矩形网格点数据、散点数据的最小平方误差双三次样条曲面拟合，e02def等可求出曲面拟合的函数值。

§5.4 范例：薄膜渗透率的测定
某种医用薄膜有允许一种物质的分子穿透它，从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜被这种分子穿透的能力。

测定方法如下：用面积为S 的薄膜将容器分成体积分别为B A V V 、的两部分，在两部分中分别注满该物质的两种不同浓度的溶液。

此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。

通过单位面积薄膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系数K 表征了薄膜被该物质分子穿透的能力，称为渗透率。

定时测量容器中薄膜某一侧的溶液浓度值，以此确定K 的数值。

5．4．1 数学模型的建立
这是一个综合性的应用实例。

主要涉及微分方程和数据拟合参数（即参数辨识）的数学知识。

1．假设：
1） 薄膜两侧的溶液始终是均匀的，即在任何时刻膜两侧的每一处溶液的浓度都是相同的。

2） 当两侧浓度不一致时，物质的分子穿透薄膜总是从高浓度溶液向低浓度溶液扩散。

3） 通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比。

4） 薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。

图4 圆柱体容器被薄膜截面S 阻隔
不同的假设应该具有不同形式的数学模型。

2．符号说明：
1）)()(t C t C B A 、表示t 时刻膜两侧溶液的浓度;
2）B A αα、表示初始时刻两侧溶液的浓度 (单位:毫克/立方厘米); 3）K 表示渗透率;
4）B A V V 、表示由薄膜阻隔的容器两侧的体积; 3．分析：
考察时段[t ，t+Δt ]薄膜两侧容器中该物质质量的变化。

以容器A 侧为例，在该时段物质质量的增加量： )()(t C V t t C V A A A A -∆+
另一方面由渗透率的定义我们知道，从B 侧渗透至A 侧的该物质的质量为：t C C SK A B ∆-)(。

由质量守恒定律，两者应该相等，于是有
)()(t C V t t C V A A A A -∆+=t C C SK A B ∆-)(。

两边除以Δt ，令Δt →0并整理得
（5.1）
且注意到整个容器的溶液中含有该物质的质量应该不变，即有下式成立：B B A A B B A A V V t C V t C V αα+=+)()(
代入（5.1）得：
再利用初始条件 ,)0(B B C α=
解出：
至此，问题归结为利用),...,2,1(N j C t C j j B =的测量数据
在时刻来辨识参数K 和B A αα,，对应的数学模型变为求函数：
∑⇒-=j
j j B B A C t C K E min ))2
（（），，（αα
令
问题转化为求函数
∑=+--+=n
j j t V V SK B A C be
a K E j B
A 1
2)1
1(
][），，（αα
的最小值点（K ，a ，b ）。

5.4.2 求解参数
例如，设1000==B A V V 立方厘米，S =10平方厘米，对容器的B 部分溶液浓度的测试结果如表1。

表1
其中C j 的单位：毫克/立方厘米。

此时极小化的函数为: ∑=⋅--+=10
1
202.0][j j t K C be a b a K E j
）
，，（
用MATLAB 软件进行计算。

1） 编写M-文件 curvefun.m
function f = curvefun(x,tdata) f = x(1)+x(2)*exp(-0.02*x(3)*tdata); 其中 x(1) = a;x(2) = b;x(3) = k; 2）编写程序(test1.m)
tdata = linspace(100,1000,10);
cdata = 1e-05.*[454 499 535 565 590 610 626 639 650 659]; x0 = [0.2,0.05,0.05];
x = curvefit('curvefun',x0,tdata,cdata)
f= curvefun(x,tdata)
e=f-cdata
plot(tdata,cdata,'o',tdata,f,'r-')
3) 输出结果:
x =
0.0070 -0.0030 0.1012
即表示k = 0.1012, a = 0.007, b = -0.003
f =
Columns 1 through 8
0.0045 0.0050 0.0054 0.0057 0.0059 0.0061 0.0063 0.0064
Columns 9 through 10
0.650.0066
误差值为：
e =
1.0e-005 *
Columns 1 through 8
-0.0339 -0.2146 0.3900 0.2857 -0.2981 -0.3570 -0.0707 0.2305
Columns 9 through 10
0.0938 -0.0339
现在curvefit函数已经被lsqcurvefit函数取代
tdata = linspace(100,1000,10);
cdata = 1e-05.*[454 499 535 565 590 610 626 639 650 659];
x0 = [0.2,0.05,0.05];
x = lsqcurvefit('curvefun',x0,tdata,cdata)
f= curvefun(x,tdata)
e=f-cdata
plot(tdata,cdata,'o',tdata,f,'r-')
运行结果
x =
0.0063 -0.0034 0.2542
f =
Columns 1 through 8
0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063
Columns 9 through 10
0.630.0063
e =
1.0e-003 *
Columns 1 through 8
-0.2322 0.1243 0.2495 0.2413 0.1668 0.0724 -0.0241 -0.1159
Columns 9 through 10
-0.2030 -0.2792
(该程
序在MATLAB 命令窗口运行)
进一步求得：004.0=B α（毫克/立方厘米），01.0=A α（毫克/立方厘米）
在MATLAB 中,还有一种求解非线性数据拟合的语句。

标准问题: min f (x ) = f 1(x )2 + f 2(x )2 + … +f m (x )2 + L 其中 L 为常数。

基本格式: x = leastsq (‘fun ’,x0) 其中fun 为f i (x )(i = 1,2,…,m)函数的M-文件。

§5.5 数据插值与拟合实验
一、实验目的及意义
[1] 了解插值、最小二乘拟合的基本原理
[2] 掌握用MATLAB计算一维插值和两种二维插值的方法；
[3] 掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法。

二、实验内容
1．针对实际问题，试建立数学模型。

用MATLAB计算一维插值和两种二维插值的方法求解；
1．用MATLAB中的函数作一元函数的多项式拟合与曲线拟合，作出误差图；
2．用MATLAB中的函数作二元函数的最小二乘拟合，作出误差图；
3．针对预测和确定参数的实际问题，建立数学模型，并求解。

三、实验步骤
1．开启软件平台——MATLAB，开启MATLAB编辑窗口；
2．根据各种数值解法步骤编写M文件
3．保存文件并运行；
4．观察运行结果(数值或图形)；
5．根据观察到的结果写出实验报告，并浅谈学习心得体会。

四、实验要求与任务
根据实验内容和步骤，完成以下具体实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论→心得体会）
1．数据插值：
（1）轮船的甲板成近似半椭圆面形
为了得到甲板的面积，首先测量得到横向最大相间8.534米；然后等间距地测得纵向高度，自左向右分别为：
0.914, 5.060, 7.772, 8.717, 9.083, 9.144, 9.083, 8.992, 8.687, 7.376,
2.073，
计算甲板的面积。

(2) 山区地貌图
在某山区（平面区域（0，2800）（0，2400）内，单位：米）测得一些地点的高程（单位：米）如下表所示，试作出该山区的地貌图和等高线图。

2．数据拟合Malthus人口指数增长模型中参数
从1790—1980年间美国每隔10年的人口记录如下表：
用以上数据检验马尔萨斯(Malthus)人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。

提示：Malthus 模型的基本假设是：人口的增长率为常数，记为 r。

记时刻t的人口为x(t），（即x(t）为模型的状态变量）且初始时刻的人口为x0，于是得到如下微分方程：
需要先求微分方程的解，再用数据拟合模型中的参数。