拟合多项式

C++代码示例如下

#include <iostream>
#include <vector>
#include<cmath>
#include <ctime>//eigen核心部分
#include <Eigen/Core>//稠密矩阵的代数运算（逆、特征值等）
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
/// <summary>
/// 拟合多项式
/// </summary>
/// <param name="degree">多项式的次</param>
/// <param name="x">x方向的数据点</param>
/// <param name="y">y方向的数据点</param>
/// <param name="a">多项式的系数</param>
/// <param name="R2">决定系数</param>
void Polynomial(int degree, std::vector<double> x, std::vector<double> y, std::vector<double> *a,double* R2) {
    // 生成矩阵A
    MatrixXd A(degree + 1, degree + 1);
    A.setZero();
    for (int j = 0; j < degree + 1; j++) {
        for (int i = 0; i < degree + 1; i++) {
            for (int k = 0; k < (int)x.size(); k++) {
                A(i, j) += 2 * pow(x[k], degree - j) * pow(x[k], degree - i);
            }
        }
    }
    // 生成向量B
    VectorXd B(degree + 1);
    B.setZero();
    for (int j = 0; j < degree + 1; j++) {
        for (int k = 0; k < (int)y.size(); k++) {
            B(j) += 2 * pow(x[k], degree - j) * y[k];
        }
    }

    //多项式方程组系数向量 = 矩阵A的逆 * 向量B
    VectorXd X = A.inverse() * B;
    for (int i = 0; i < X.size(); i++) {
        a->push_back(X(i));//把多项式的系数推到向量指针里
    }

    //计算决定系数R2
    double sum_y = 0.0;
    for (int i = 0; i < (int)y.size(); i++) {
        sum_y += y[i];
    }
    double mean_y = sum_y / y.size();//计算平均值
    double ssr = 0, sst = 0;
    std::vector<double> fit_y;//拟合y
    for (int i = 0; i < (int)y.size(); i++) {
        double value = 0;
        for (int n = degree; n >= 0; n--) {
            value += X(degree - n) * pow(x[i], n);//逐个计算拟合后的y
        }
        fit_y.push_back(value);//把拟合后的y推到向量里
    }
    for (int i = 0; i < (int)y.size(); i++) {
        ssr += pow(fit_y[i] - mean_y, 2);
        sst += pow(y[i] - mean_y, 2);
    }
    *R2 = ssr / sst;
}
int main(int argc, char** argv) {
    // 定义多项式的次
    int degree = 3;

    // 定义数据点
    std::vector<double> x = { 1, 2, 3 ,4, 5 };
    std::vector<double> y = { 1.1, 3.2, 6.7, 18.3, 100.3 };
    std::vector<double>* a = new std::vector<double>();;
    double R2 = 0;
    double* ptr = &R2;
    Polynomial(degree,x,y,a, ptr);
    //打印多项式系数 a0 a1... an
    int i = 0;
    for (double elem : *a) {
        cout << "a" << i++ << " = " << elem << endl;
    }
    cout << "R2 = "<< R2 << endl;
}

main.cpp

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

namespace QuadraticFitExample
{
    class Program
    {
        static void Main(string[] args)
        {
            // 假设我们有一些样本数据点
            double[] xData = { 1, 2, 3, 4, 5 };
            double[] yData = { 1.1, 3.2, 6.7, 18.3, 100.3 };

            // 指定多项式的阶数（这里是3次多项式）
            int degree = 3;

            // 使用LeastSquaresFit进行多项式拟合
            double[] coefficients = Fit.Polynomial(xData, yData, degree);

            Console.WriteLine("多项式系数:");
            for (int i = 0; i <= degree; i++)
            {
                Console.WriteLine($"x^{i}项系数: {coefficients[i]}");
            }

            // 计算预测值（考虑到所有五次项）
            double[] predictedValues = xData.Select(x =>
                coefficients.Select((coeff, power) => coeff * Math.Pow(x, power)).Sum()).ToArray();

            // 计算均方误差和总平方和
            double mse = MeanSquaredError(yData, predictedValues);
            double tss = TotalSumOfSquares(yData, Enumerable.Average(yData));

            // 计算决定系数 R^2
            double rSquared = 1.0 - mse / tss;

            Console.WriteLine($"\n决定系数 R^2: {rSquared}");

            // 验证拟合效果（对原数据点进行预测并打印）
            foreach (var x in xData.Zip(predictedValues, (x, pred) => (x, pred)))
            {
                Console.WriteLine($"Predicted y for x={x.x}: {x.pred}");
            }

            Console.ReadKey();
        }

        // 辅助方法计算均方误差
        static double MeanSquaredError(double[] actuals, double[] predictions)
        {
            return actuals.Zip(predictions, (a, p) => Math.Pow(a - p, 2)).Average();
        }

        // 辅助方法计算总平方和
        static double TotalSumOfSquares(double[] values, double mean)
        {
            return values.Sum(v => Math.Pow(v - mean, 2));
        }
    }
}