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

线性回归

让我们首先考虑最小二乘拟合的线性形式

其中fi(x)是已经定义过的x的函数，叫做基本函数。这种情况下，残差和就是

相应的梯度方程从简化到

先前的方程能写成矩阵形式

这个矩阵形式也叫做最小二乘拟合的标准方程，它可以通过前面章节讨论的Gauss-Jordan方法求解。

 

实现

使用前面章节提到的算法，我们可以实现线性回归方法。添加一个新的公有静态方法LinearRegression到CurveFittingAlgorithms类里：

public delegate doubleModelFunction(double x);public static VectorRLinearRegression(double[] xarray, double[] yarray,                                      ModelFunction[] f, out double sigma){    int m = f.Length;    MatrixR A = new MatrixR(m, m);    VectorR b = new VectorR(m);    int n = xarray.Length;     for (int k = 0; k < m; k++)    {        b[k] = 0.0;        for (int i = 0; i < n; i++)        {            b[k] += f[k](xarray[i]) *yarray[i];        }    }     for (int j = 0; j < m; j++)    {        for (int k = 0; k < m; k++)        {            A[j, k] = 0.0;            for (int i = 0; i < n; i++)            {                A[j, k] += f[j](xarray[i]) *f[k](xarray[i]);            }        }    }    VectorR coef = GaussJordan(A, b);     // Calculate the standard deviation:    double s = 0.0;    for (int i = 0; i < n; i++)    {        double s1 = 0.0;        for (int j = 0; j < m; j++)        {            s1 += coef[j] * f[j](xarray[i]);        }        s += (yarray[i] - s1) * (yarray[i] -s1);    }    sigma = Math.Sqrt(s / (n - m));    return coef;}

注意到我们开始定义了一个委托函数，它将一个双精度变量x作为输入参数。接着我们实现了一个公有精通方法LinearRegression，它返回一个VectorR对象并将其作为基本方程的系数。

       在LinearRegression方法里，我们指定系数矩阵A和向量b。我们接着通过GaussJordan方法求解标准方程.根据解决方案,我们计算了标准差.这个方法返回多项式(一个VectorR的对象)的系数;我们在PolynomialFit方法里使用out前缀指定标准差。

public static VectorRLinearRegression(double[] xarray, double[] yarray,                      ModelFunction[] f, outdouble sigma){     }

这样，你可以从一个C#方法获得多个输出结果。

 

测试线性回归

你可以简单地使用前面章节提到的方法展示多项式拟合。添加一个新的WPF Windows 应用程序到当前项目，命名为PolynomialFit.下面是这个例子的XAML文件：

<Windowx:Class="CurveFitting.PolynomialFit"   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="PolynomialFit"Height="450" Width="500">    <Grid x:Name="rootGrid"SizeChanged="rootGrid_SizeChanged">        <chart:LineChartControlLibx:Name="myChart" Xmin="0" Xmax="6"XTick="1"               Ymin="-100"Ymax="600" YTick="100" Title="PolynomialFitting"/>    </Grid></Window>

这个XAML文件创建了一个叫做myChart的线图表。我们将使用它来展示曲线拟合的结果。下面是相关的后台代码：

using System;usingSystem.Windows;usingSystem.Windows.Controls;usingSystem.Windows.Media; namespaceCurveFitting{    public partial class PolynomialFit : Window    {        public PolynomialFit()        {            InitializeComponent();        }        private voidrootGrid_SizeChanged(object sender, SizeChangedEventArgs e)        {            myChart.Width =rootGrid.ActualWidth;            myChart.Height =rootGrid.ActualHeight;            AddData();        }         private void AddData()        {            double[] x0 = new double[] { 1, 2,3, 4, 5 };            double[] y0 = new double[] { 5.5,43.1, 128, 290.7, 498.4 };             VectorR[] results = new VectorR[3];             for (int i = 0; i <results.Length; i++)            {                double sigma = 0;                results[i] =CurveFittingAlgorithms.PolynomialFit(x0, y0,                             i + 1, out sigma);            }             // Plot results:           myChart.DataCollection.DataList.Clear();            LineCharts.DataSeries ds;           // Plot results:           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.Triangle;           ds.Symbols.BorderColor =Brushes.Black;           for (int i = 0; i < x0.Length;i++)           {               ds.LineSeries.Points.Add(newPoint(x0[i], y0[i]));           }           myChart.DataCollection.DataList.Add(ds);           // 1st order fitting data:           ds = new LineCharts.DataSeries();           ds.LineColor = Brushes.DarkGreen;           ds.LineThickness = 2;           ds.SeriesName = "1st OrderFitting";           for (int i = 0; i < 141; i++)           {               double x = -1.0 + i / 20.0;               double y = results[0][0] +results[0][1] * x;               ds.LineSeries.Points.Add(newPoint(x, y));           }          myChart.DataCollection.DataList.Add(ds);           // 2nd order fitting data:           ds = new LineCharts.DataSeries();           ds.LineColor = Brushes.Red;            ds.LineThickness = 2;            ds.LinePattern = LineCharts.DataSeries.LinePatternEnum.Dash;            ds.SeriesName = "2nd OrderFitting";            for (int i = 0; i < 141; i++)            {                double x = -1.0 + i / 20.0;                double y = results[1][0] +results[1][1] * x +                           results[1][2] * x * x;                ds.LineSeries.Points.Add(newPoint(x, y));            }           myChart.DataCollection.DataList.Add(ds);             // 3rd order fitting data:            ds = new LineCharts.DataSeries();            ds.LineColor = Brushes.DarkBlue;            ds.LineThickness = 2;            ds.LinePattern =LineCharts.DataSeries.LinePatternEnum.DashDot;            ds.SeriesName = "3rd OrderFitting";            for (int i = 0; i < 141; i++)           {                double x = -1.0 + i / 20.0;                double y = results[2][0] +results[2][1] * x +                           results[2][2] * x *x + results[2][3] * x * x * x;                ds.LineSeries.Points.Add(newPoint(x, y));            }           myChart.DataCollection.DataList.Add(ds);            myChart.IsLegend = true;            myChart.LegendPosition =LineCharts.Legend.LegendPositionEnum.NorthWest;        }    }}

这里，我们对数据点拟合了一个不同阶数（从1阶到3阶）的的多项式。运行程序生成的结果如图11-3所示。你可以从图中看到对这个指定的数据集，二阶多项式的拟合结果与三阶多项式的拟合结果几乎是相同的。

图11-3多项式曲线拟合

 

加权线性回归

我们前面讨论到，线性回归对于通过模型函数展示观察数据是一种有用的技术，可以表示成最小二乘法的最小化问题。在线性最小二乘中，结果分析需要使用Gauss-Jordan方法求解联立方程组。线性回归中的一个重要的假设是所有的错误具有相同的意义。

然而，有些情况下不同的点的精确度的可靠性不一样。例如，有可能存在一个测量的精度漂移，一些错误或多或少比其他更重要。为了考虑到这些因素，我们引入加权因子到每个数据点来最小化加权的残差平方和：

这个过程强制将拟合函数f(x)朝更大权重的点逼近。

       对于最简单的线性回归，例如，直线拟合，它的拟合函数通过f(x)=a+bx给出，前面的方程就变成

根据最小化约束条件和，我们可以确定系数a和b：

这里，xw和yw是加权平均：

 

实现

使用前面提出的方法，我们可以实现加权线性回归。添加一个新的公有静态方法WeightedLinearRegression到CurveFittingAlgorithms类

public static double[]WeightedLinearRegression(double[] xarray, double[] yarray, double[] warray){    int n = xarray.Length;    double xw = 0.0;    double yw = 0.0;    double b1 = 0.0;    double b2 = 0.0;    double a = 0.0;    double b = 0.0;     for (int i = 0; i < n; i++)    {        xw += xarray[i] / n;        yw += yarray[i] / n;    }    for (int i = 0; i < n; i++)    {        b1 += warray[i] * warray[i] * yarray[i]* (xarray[i] - xw);        b2 += warray[i] * warray[i] * xarray[i]* (xarray[i] - xw);    }    b = b1 / b2;    a = yw - xw * b;    return new double[] { a, b };}

这里xarray和yarray是x和y数据的集合，表示给定的数据点。Warray是每个数据点的权重。这个方法返回模型函数的系数a和b：

 

指数函数拟合

我们可以使用加权线性回归对数据点进行指数函数拟合。例如，对模型函数，当在最小二乘技术用使用时，通常会导致

系数a和b的非线性依赖。然而，如果我们使用logy而不是y，问题就变成了线性回归。这种情况下，拟合函数变成

       注意到数据对数的最小二乘拟合与原始数据的最小二乘拟合不太一样。对数拟合的残差是

而原始数据拟合的残差是

我们可以对对数拟合加权来消除这种差异。根据前面的方程，我们有

因此对数拟合的残差能够重新写成下面的形式：

在约束中，我们可以使用近似处理。你可以通过最小化来明白这个，我们需要引入权重因子。当对数据点拟合F(x)时如果我们使用权重我们可以这样做。因此，最小化

将变成最小化的近似。

       加上我们有下面的一组数据：

x = 1, 2 , 3, 4, 5 , 6, 7,8, 9, 10

y = 1.9398, 2.9836,5.9890, 10.2, 20.7414, 23.232, 69.5855, 82.5836, 98.1779 ,339.3256

我们想要使用指数函数来拟合这些数据。首先我们需要进行一个对数变换log y = log a +bx，然后使用WeightedLinearRegression方法计算loga和b的系数。

       添加一个新的WPF Windows应用程序到当前项目中，命名为WeightedLinearRegression.下面是这个例子的XAML文件：

<Windowx:Class="CurveFitting.WeightedLinearRegression" 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="WeightedLinearRegression"Height="320" Width="600">    <Grid x:Name="rootGrid"SizeChanged="rootGrid_SizeChanged">        <StackPanelOrientation="Horizontal">            <chart:LineChartControlLibx:Name="logChart"                   Xmin="0"Xmax="11" XTick="1" Title="Log Scale Plot"/>            <chart:LineChartControlLibx:Name="linearChart"                   Xmin="0"Xmax="11" XTick="1" Title="Linear ScalePlot"/>        </StackPanel>    </Grid></Window>

这里我们创建了两个线图表控件:一个用来显示log缩放的结果，另一个用来显示线性缩放的结果。下面是相关的后台代码：

using System;using System.Windows;usingSystem.Windows.Controls;usingSystem.Windows.Media; namespace CurveFitting{    public partial classWeightedLinearRegression : Window    {        public WeightedLinearRegression()        {            InitializeComponent();        }         private voidrootGrid_SizeChanged(object sender, SizeChangedEventArgs e)        {            logChart.Width =rootGrid.ActualWidth / 2;            logChart.Height =rootGrid.ActualHeight;            linearChart.Width =rootGrid.ActualWidth / 2;            linearChart.Height =rootGrid.ActualHeight;            AddData();        }     private void AddData()        {            double[] x0 = new double[] { 1, 2,3, 4, 5, 6, 7, 8, 9, 10 };            double[] y0 = new double[] {1.9398, 2.9836, 5.9890, 10.2, 20.7414,                23.232, 69.5855, 82.5836,98.1779, 339.3256 };            double[] ylog = new double[] {0.6626, 1.0931, 1.7899, 2.3224,                3.0321, 3.1455, 4.2426, 4.4138,4.5868, 5.8270 };             double[] results =               CurveFittingAlgorithms.WeightedLinearRegression(x0, ylog, y0);             // Plot linear scale results:            LinearScale(x0, y0, results);             // Plot log scale results:         LogScale(x0, ylog, results);     }     private void LinearScale(double[] x0,double[] y0, double[] results)     {        linearChart.DataCollection.DataList.Clear();         LineCharts.DataSeries ds;         linearChart.ChartStyle.Ymin = -50;         linearChart.ChartStyle.Ymax = 350;         linearChart.ChartStyle.YTick = 50;         linearChart.ChartStyle.YLabel ="Y";         // Plot original data:         ds = new LineCharts.DataSeries();         ds.LineColor = Brushes.Transparent;         ds.SeriesName = "Original";         ds.Symbols.SymbolType =LineCharts.Symbols.SymbolTypeEnum.Triangle;         ds.Symbols.BorderColor =Brushes.Black;         for (int i = 0; i < x0.Length; i++)         {             ds.LineSeries.Points.Add(newPoint(x0[i], y0[i]));         }        linearChart.DataCollection.DataList.Add(ds);         // Plot curve fittin data:         ds = new LineCharts.DataSeries();         ds.LineColor = Brushes.DarkGreen;         ds.LineThickness = 2;         ds.SeriesName = "CurveFitting";         for (int i = 0; i < 111; i++)         {             double x = 0.1 + i / 10.0;             double y = Math.Exp(results[0] +results[1] * x);             ds.LineSeries.Points.Add(newPoint(x, y));         }        linearChart.DataCollection.DataList.Add(ds);         linearChart.Legend.IsLegend = true;         linearChart.Legend.LegendPosition =            LineCharts.Legend.LegendPositionEnum.NorthWest;     }   private void LogScale(double[] x0, double[]ylog, double[] results)   {      logChart.DataCollection.DataList.Clear();       LineCharts.DataSeries ds;       logChart.ChartStyle.Ymin = 0;       logChart.ChartStyle.Ymax = 7;       logChart.ChartStyle.YTick = 1;       logChart.ChartStyle.YLabel ="Log(y)";       // Plot original data:            ds = new LineCharts.DataSeries();            ds.LineColor = Brushes.Transparent;            ds.SeriesName ="Original";            ds.Symbols.SymbolType = LineCharts.Symbols.SymbolTypeEnum.Triangle;            ds.Symbols.BorderColor =Brushes.Black;            for (int i = 0; i < x0.Length;i++)            {                ds.LineSeries.Points.Add(newPoint(x0[i], ylog[i]));            }            logChart.DataCollection.DataList.Add(ds);             // Plot curve fittin data:            ds = new LineCharts.DataSeries();            ds.LineColor = Brushes.DarkGreen;            ds.LineThickness = 2;            ds.SeriesName = "CurveFitting";            for (int i = 0; i < 111; i++)            {                double x = 0.1 + i / 10.0;                double y = results[0] +results[1] * x;                ds.LineSeries.Points.Add(newPoint(x, y));            }            logChart.DataCollection.DataList.Add(ds);            logChart.Legend.IsLegend = true;            logChart.Legend.LegendPosition =               LineCharts.Legend.LegendPositionEnum.NorthWest;        }    }}

       在AddData方法里，语句

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

将生成结果( -0.0686983921044781,0.578232434928087)。第一个数字是(loga)的参数，而第二个是参数b，他们得到

       我们可以通过将它们绘制出来来检验拟合结果的好坏（如图11-4）.你可以从图11-4中看到加权线性回归确实给出了一个很好的拟合结果。

图11-4加权线性回归的结果：对数绘图(左)和线性绘图(右)