《Practical WPF Charts and Graphics 》翻译——之11章 曲线拟合(1)
第11章
曲线拟合
在科学和工程里,实验得到的数据通常因为测量误差会包含许多随机噪声。曲线拟合的目的是找到一个光滑的曲线去拟合平均数据点。我们通常要求这个曲线具有简单低阶多项式的形式,才能不再产生数据的随机误差。
插值和拟合有一定的区别。插值,在前一章讨论过,可以看着是曲线拟合的特殊形式,它的函数必须精确穿过数据点。这就暗示给定的插值数据是精确和明显的。曲线拟合应用于通常因为实验测量误差而包含噪声的数据。它试图找到给定数据集的最佳拟合。因此,曲线不需要穿过给定的数据点。
线性代数方程
曲线拟合技术经常需要解线性代数方程。在本章,我们会回顾一种常用的方法,Gauss-Jordan消元法,用来解线性方程。通常,线性方程系统包含多个变量。广义向量和矩阵分析是线性系统的基础。
注意WPF包含的向量和矩阵结构被设计专门用来对各种图形元素进行操作的,但是他们对于线性方程系统这样的应用程序不够通用。例如,如果一个线性系统包含10个方程和10个未知数,系统里的向量和矩阵将会是10维的。 WPF里的2DVector和Matrix或者3D Vector3D和Matrix3D结构不能应用于这样的系统。因此,我们需要创建通用n维向量和n*n矩阵结构。
在我的另一本书里,Practical Numerical Methods with C#,我详细展示了怎么创建这样的用于实变量和复变量的广义向量和矩阵。这里,我不推导而直接使用VectorR和MatrixR结构用于实变量。如果你对这些结构的实现感兴趣,请参考Practical Numerical Methods with C#。
这里,我们将考虑怎么求解包含任意数量的未知数的线性方程。一系列线性方程能够写成下面的矩阵形式
这里,A是矩阵系数的方矩阵,x和b都是表示未知数和右边常数的列向量。在下面的部分里,我将展示使用Gauss-Jordan方法解线性方程。
Gauss-Jordan算法
使用Gauss-Jordan消元法求解线性方程组会得到方程的解和系数矩阵的逆。Gauss-Jordan消元法的过程只需要两步。第一步,叫做向前消元,将一个给定的系统简化成三角或者梯形形式或者退化的结果来表示系统无解。我们通过使用基本操作来完成它。第二步,叫做向后消元,使用回代来找到线性方程的解。
就矩阵形式而言,第一步使用基本操作将一个矩阵简化到行梯阵的形式,第二部将矩阵减少到行标准形式。
Gauss-Jordan消元法使用基本操作计算矩阵分解:行相乘,行交换,和向其他行加上某一行的倍数。算法的第一部分计算分解;第二部分将原始矩阵写成唯一确定的可逆矩阵和唯一确定简化行梯阵式的乘积。Gauss-Jordan方法通常和其他技术效率一样而且十分稳定。
Gauss-Jordan消元法最简单的方法是从所有i>0里的方程中消去方程(11.1)里的x[0],从所有i>1的方程中消去x[1];一般来说,从i>j的方程里消去x[j]。这个过程,叫做Gauss-Jordan消元,生成主对角线下方系数全为零的三角矩阵。
使用这个方法我们必须注意。要是在三角化得过程中一个对角元素变成零或者太小,这个方法就会失败。我们可以通过互换这两套方程来避免这个问题。方程表示的顺序并不影响它们的解,但是它对你怎么求解这个方程组会有很大影响。因此,我们重新排列方程组来使对角元素尽可能大。我们会使用Pivot方法来完成这个任务。
实现
新建一个WPF Windows应用程序,命名为CurveFitting.首先我们需要通过鼠标右键单击,选择Add|Add Existing Item…向当前项目中添加VectorR和MatrixR结构。你可以打开VectorR.cs和MatrixR.cs文件来查看源代码,看他们怎么实现的。你可以在Practical Numerical Methods with C#里找到它们实现的详细解释。给项目添加一个新类,命名为CurveFittingAlgorithms.这里我们只考虑线性系统和实数变量。下面是类的代码。
这个函数是要使用有n+1个数据点的数据集来拟合的。前面的模型函数包含m+1个变量参数a0,a1,….,am,其中m<n。我们希望找到模型最佳拟合数据的参数值。
模型函数的形式在前面定义了,通常是从获取观察数据的实验的相关理论中来的。曲线拟合过程包含两个步骤:选择模型函数的形式,接下来是计算数据最佳拟合的参数。
最小二乘拟合使残差平方和最小
相对于每个参数ai。方程的术语叫做残差,它表示拟合函数和数据点的差异。因此,参数的最优值由下面的条件给定:
前面的方程对ai来说一般是非线性的而且很难求解。这种梯度方程适用于所有的最小二乘问题。每一个特别的问题需要模型函数的特别的表达式和它的偏导数。例如,模型函数经常选择为特定函数fi(x)的线性组合:
这种情况下,梯度方程就变成线性的了。如果模型函数是一个多项式,我们有f0(x) = 1, f1(x) = x, f2(x) = x2等等。
直线拟合
最简单的线性回归是直线拟合,它试图使用最小二乘拟合一条直线。它考察一个独立的变量x和一个独立的变量y之间的相关性。这种情况下,模型函数有下面的简单形式:
线性回归要最小化的和函数就变成
相应的梯度方程变为
我们能够通过前面的两个方程找到参数a和b的一个解:
这里,xm和ym是x和y数据的均值:
标准差σ能表示为
实现
添加一个新的方法到CurveFittingAlgorithms类,命名为StraightLineFit。下面是方法的代码:
而且标准差是0.191。
图11-1显示了原始数据直线拟合的结果。
图11-1 直线拟合
(接下文——线性回归)