C#实战：如何用最小二乘法将两组数据拟合成二次函数

1次阅读

共计 2113 个字符，预计需要花费 6 分钟才能阅读完成。

在工程预测和实验分析中，我们常常需要根据观测数据建立数学模型。比如根据历史销售数据预测未来趋势，或者通过实验测量值推导物理规律。数据拟合就是将离散的观测点转化为连续函数的过程，而二次函数因其良好的平滑性和计算简便性，成为常用的拟合模型之一。

拟合二次函数主要有两种方法：解析法和数值法。解析法直接求解方程组，但当数据存在噪声或矩阵条件数较大时容易失效。相比之下，最小二乘法这种数值方法通过最小化误差平方和来求解，具有更好的数值稳定性。其数学表达为：

$$
\min \sum_{i=1}^n (y_i – (ax_i^2 + bx_i + c))^2
$$

首先需要安装 Math.NET Numerics 包。在 NuGet 包管理器控制台中运行：

Install-Package MathNet.Numerics

建立范德蒙德矩阵的关键代码：

// 数据标准化预处理
double[] xValues = {1, 2, 3, 4, 5};
double[] yValues = {1.2, 3.9, 8.1, 15.3, 24.5};

// 构建范德蒙德矩阵
var A = Matrix<double>.Build.DenseOfRowArrays(xValues.Select(x => new[] {x*x, x, 1}).ToArray());
var b = Vector<double>.Build.Dense(yValues);

当数据点分布不均匀时，矩阵可能呈现病态。这时可以引入 Tikhonov 正则化：

// 正则化参数
double lambda = 0.1;
var I = Matrix<double>.Build.DenseIdentity(3);

// 正则化后的解
var coefficients = (A.Transpose() * A + lambda * I).Inverse() * A.Transpose() * b;

评估拟合质量的重要指标：

double yMean = yValues.Average();
double ssTot = yValues.Sum(y => Math.Pow(y - yMean, 2));
double ssRes = yValues.Zip(A * coefficients, (y, yFit) => Math.Pow(y - yFit, 2)).Sum();
double rSquared = 1 - ssRes/ssTot;

数据量不足时的处理 ：
当数据点少于 3 个时，建议增加数据或降低多项式阶数
可以尝试添加虚拟约束点保持矩阵非奇异
浮点数精度问题 ：
对数据进行标准化处理
使用 decimal 类型进行关键计算
避免直接比较浮点数，改用容差比较
可视化验证 ：
使用 ScottPlot 或 OxyPlot 绘制原始数据和拟合曲线
观察残差分布是否随机

using MathNet.Numerics.LinearAlgebra;
using System;
using System.Linq;

class Program
{static void Main()
    {
        try
        {
            // 输入数据
            double[] xData = {1, 2, 3, 4, 5};
            double[] yData = {1.2, 3.9, 8.1, 15.3, 24.5};

            // 构建矩阵
            var A = Matrix<double>.Build.DenseOfRowArrays(xData.Select(x => new[] {x*x, x, 1}).ToArray());
            var b = Vector<double>.Build.Dense(yData);

            // 解最小二乘问题
            var coefficients = A.TransposeThisAndMultiply(A).Inverse() * A.Transpose() * b;

            Console.WriteLine($"拟合方程: y = {coefficients[0]:F2}x² + {coefficients[1]:F2}x + {coefficients[2]:F2}");

            // 计算 R²
            double yMean = yData.Average();
            double ssTot = yData.Sum(y => Math.Pow(y - yMean, 2));
            double ssRes = yData.Zip(A * coefficients, (y, yFit) => Math.Pow(y - yFit, 2)).Sum();
            Console.WriteLine($"R² = {1 - ssRes/ssTot:F4}");
        }
        catch(SingularMatrixException)
        {Console.WriteLine("矩阵奇异，请检查数据或添加正则化项");
        }
    }
}