《Practical WPF Charts and Graphics 》翻译——之11章 曲线拟合（1）

第11章

曲线拟合

 

在科学和工程里，实验得到的数据通常因为测量误差会包含许多随机噪声。曲线拟合的目的是找到一个光滑的曲线去拟合平均数据点。我们通常要求这个曲线具有简单低阶多项式的形式，才能不再产生数据的随机误差。

插值和拟合有一定的区别。插值，在前一章讨论过，可以看着是曲线拟合的特殊形式，它的函数必须精确穿过数据点。这就暗示给定的插值数据是精确和明显的。曲线拟合应用于通常因为实验测量误差而包含噪声的数据。它试图找到给定数据集的最佳拟合。因此，曲线不需要穿过给定的数据点。

 

线性代数方程

 

曲线拟合技术经常需要解线性代数方程。在本章，我们会回顾一种常用的方法，Gauss-Jordan消元法，用来解线性方程。通常，线性方程系统包含多个变量。广义向量和矩阵分析是线性系统的基础。

注意WPF包含的向量和矩阵结构被设计专门用来对各种图形元素进行操作的，但是他们对于线性方程系统这样的应用程序不够通用。例如，如果一个线性系统包含10个方程和10个未知数，系统里的向量和矩阵将会是10维的。 WPF里的2D Vector和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.这里我们只考虑线性系统和实数变量。下面是类的代码。

using System;
using System.Collections;
namespace CurveFitting
{
    public classCurveFittingAlgorithms
    {
        public staticVectorR GaussJordan(MatrixR A, VectorR b)
        {
            Triangulate(A,b);
            int n =b.GetSize();
            VectorR x =new VectorR(n);
            for (int i = n- 1; i >= 0; i--)
            {
                double d =A[i, i];
                if(Math.Abs(d) < 1.0e-500)
                    thrownew ArgumentException("Diagonal element is too small!");
                x[i] =(b[i] - VectorR.DotProduct(A.GetRowVector(i), x)) / d;
          }
          return x;
      }
      private static voidTriangulate(MatrixR A, VectorR b)
      {
          int n =A.GetRows();
          VectorR v = newVectorR(n);
          for (int i = 0;i < n - 1; i++)
          {
              double d =Pivot(A, b, i);
              if(Math.Abs(d) < 1.0e-500)
                  thrownew ArgumentException("Diagonal element is too small!");
              for (int j =i + 1; j < n; j++)
              {
                  double dd = A[j, i] / d;
                  for (intk = i + 1; k < n; k++)
                  {
                      A[j,k] -= dd * A[i, k];
                  }
                  b[j] -=dd * b[i];
              }
          }
      }
 
        private static double Pivot(MatrixR A,VectorR b, int q)
        {
            int n =b.GetSize();
            int i = q;
            double d =0.0;
            for (int j =q; j < n; j++)
            {
                double dd= Math.Abs(A[j, q]);
                if (dd> d)
                {
                    d =dd;
                    i = j;
                }
            }
            if (i > q)
            {
               A.GetRowSwap(q, i);
               b.GetSwap(q, i);
            }
            return A[q,q];
        }
    }
}

       这个类包含实现Gauss-Jordan消元法的三个方法。Pivot方法通过重排方程来返回最大可用对角元素。如果需要，我们通过使用在MatrixR和VectorR里定义的Swap方法交换系数矩阵A和常向量b的行。

       我们通过使用Triangulate方法实现三角化。在第一层循环中，我们使用Pivot的最大可能的元素而不是直接使用对角元素A[i,i]。

       最后，我们实现了公有的GaussJordan方法。在这个方法里，我们使用回代找到了三角方程的未知数的值。这个过程从最后一个方程开始，它只包含一个未知数，意味着我们可以立即找出它的解。我们可以将这个值代入到倒数第二个方程里面来计算下一个未知数。我们重复这个过程，直到求出所有的未知数。

 

最小二乘拟合

 

最流行的曲线拟合方法是最小二乘拟合，它经常用来求解超定系统。它经常在统计学中使用，特别是回归分析。最小二乘意义上的最佳拟合是残差平方和最小的模型实例。最小二乘里的残差是指观察数据里的值和模型提供的值的差别。

       假设一组给定的数据由n+1个数据点(xi,yi),i=0,1,….n组成，其中xi是一个独立的变量，yi也是一个独立的变量，它们的值通过观测获取。模型函数定义如下：

      

这个函数是要使用有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。下面是方法的代码：

public static double[] StraightLineFit(double[] xarray, double[] yarray)
{
    int n = xarray.Length;
    double xm = 0.0;
    double ym = 0.0;
    double b1 = 0.0;
    double b2 = 0.0;
    double a = 0.0;
    double b = 0.0;
    double s = 0.0;
    double sigma = 0.0;
    for (int i = 0; i < n; i++)
    {
        xm += xarray[i] / n;
        ym += yarray[i] / n;
    }
    for (int i = 0; i < n; i++)
    {
        b1 += yarray[i] * (xarray[i] - xm);
        b2 += xarray[i] * (xarray[i] - xm);
    }
    b = b1 / b2;
    a = ym - xm * b;
    for (int i = 0; i < n; i++)
    {
        s += (yarray[i] - a - b * xarray[i]) *(yarray[i] - a - b * xarray[i]);
    }
    sigma = Math.Sqrt(s / (n - 2));
    return new double[] {a, b, sigma };
}

测试直线拟合

你可以简单地使用前面章节提到的方法来展示直线拟合。为了直接将结果显示在你的屏幕上，你需要添加LineChartControl到项目的引用，就像我们在前面章节中做的。添加一个新的WPF Windows应用程序到当前项目并命名为StraightLineFit。下面是例子的XAML文件：

<Windowx:Class="CurveFitting.StraightLineFit"
   xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
   xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
   xmlns:chart="clr-namespace:LineChartControl;assembly=LineChartControl"
    Title="StraightLineFit"Height="350" Width="400">
    <Grid x:Name="rootGrid"SizeChanged="rootGrid_SizeChanged">
        <chart:LineChartControlLibx:Name="myChart" Xmin="-1" Xmax="6"
               XTick="1"Ymin="1" Ymax="7" YTick="1" Title="StraightLine Fitting"/>
    </Grid>
</Window>

这个XAML创建了一个叫做myChart的线图表控件，它用来显示曲线拟合的结果。相关的后台代码如下：

using System;
using System.Windows;
usingSystem.Windows.Controls;
usingSystem.Windows.Media;
namespace CurveFitting
{
    public partial class StraightLineFit :Window
    {
        public StraightLineFit()
        {
            InitializeComponent();
        }
 
        private voidrootGrid_SizeChanged(object sender, SizeChangedEventArgs e)
        {
            myChart.Height =rootGrid.ActualHeight;
            myChart.Width =rootGrid.ActualWidth;
            AddData();
        }
 
        private void AddData()
        {
            double[] x0 = new double[] { 0, 1,2, 3, 4, 5 };
            double[] y0 = new double[] { 1.9, 2.7, 3.3,4.4, 5.5, 6.5 };
            double[] results =CurveFittingAlgorithms.StraightLineFit(x0, y0);
 
           myChart.DataCollection.DataList.Clear();
            LineCharts.DataSeries ds;
            
            // Plot original data
            ds = new LineCharts.DataSeries();
            ds.LineColor = Brushes.Transparent;
            ds.SeriesName ="Original";
            ds.Symbols.SymbolType=LineCharts.Symbols.SymbolTypeEnum.OpenTriangle;
            ds.Symbols.BorderColor =Brushes.DarkBlue;
            for (int i = 0; i < x0.Length;i++)
            {
                ds.LineSeries.Points.Add(newPoint(x0[i], y0[i]));
            }
           myChart.DataCollection.DataList.Add(ds);
 
            // Curve fitting data:
            ds = new LineCharts.DataSeries();
            ds.LineColor = Brushes.DarkGreen;
            ds.SeriesName = "CurveFitting";
            for (int i = 0; i < 101; i++)
            {
                double x = i / 20.0;
                double y = results[0] +results[1] * x;
                ds.LineSeries.Points.Add(newPoint(x, y));
            }
           myChart.DataCollection.DataList.Add(ds);
            myChart.IsLegend = true;
            myChart.LegendPosition =LineCharts.Legend.LegendPositionEnum.NorthWest;
        }
    }
}

这里，我们使用两个双精度数组xarray和yarray提供了输入数据，并且使用下面的语句来计算直线拟合的系数。

double[] results =CurveFittingAlgorithms.StraightLineFit(x0, y0);

       你可以通过添加下面的代码段来使用一个Message Box来查看AddData方法的结果：

MessageBox.Show(results[0].ToString()+ ", " +
                results[1].ToString() + "," +
                results[2].ToString());

它应该显示(1.729, 0.929, 0.191)。因此，这条直线能表示成

而且标准差是0.191。

图11-1显示了原始数据直线拟合的结果。

图11-1 直线拟合

 

 （接下文——线性回归）