纯C++做多项式拟合

一、多项式拟合用途

当前有一组对应的x、y数据，希望通过这些数据点做出近似的多项式曲线：Y=···+nX^2+mX+c
其中多项式最高次数可调，返回各个参数及曲线的拟合度R^2

二、函数实现

参数中的order为设置的多项式最高次次数，coefficients为各次的系数

double polynomialFit(vector<double>& x, vector<double>& y, unsigned char order, vector<double>& coefficients){
    if (x.size() <= order || y.size() <= order) {
    return 0;
    }

    // 构建矩阵A和向量b
    int m = x.size();
    int n = order + 1;
    vector<vector<double>> A(n, vector<double>(n, 0));
    vector<double> b(n, 0);

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < m; k++) {
                A[i][j] += pow(x[k], i+j);
            }
        }

        for (int k = 0; k < m; k++) {
            b[i] += y[k] * pow(x[k], i);
        }
    }

    n = A.size();

    for (int i = 0; i < n; i++) {
        // 列主元素消去
        int maxRow = i;
        double maxVal = fabs(A[i][i]);
        for (int k = i + 1; k < n; k++) {
            if (fabs(A[k][i]) > maxVal) {
                maxVal = fabs(A[k][i]);
                maxRow = k;
            }
        }
        if (maxRow != i) {
            std::swap(A[i], A[maxRow]);
            std::swap(b[i], b[maxRow]);
        }

        // 消元过程
        for (int k = i + 1; k < n; k++) {
            double factor = A[k][i] / A[i][i];
            for (int j = i; j < n; j++) {
                A[k][j] -= factor * A[i][j];
            }
            b[k] -= factor * b[i];
        }
    }

    // 回代求解
    vector<double> result(n, 0);
    for (int i = n - 1; i >= 0; i--) {
        double temp = b[i];
        for (int j = i + 1; j < n; j++) {
            temp -= A[i][j] * result[j];
        }
        result[i] = temp / A[i][i];
    }
    coefficients = result;
    double SSR=0;
    double SST=0;
    double sumY = std::accumulate(std::begin(y),std::end(y),0.0);
    double avgY = sumY/y.size();
    for(uint16_t i=0;i<y.size();i++){
        double actY=0;
        for(unsigned char j=0;j<=order;j++){
            actY+=pow(x[i],j)*result[j];
        }
        SSR += pow(actY-y[i],2);
        SST += pow(y[i]-avgY,2);
    }
    double R = 1-SSR/SST;
    return R;
}

三、函数调用

vector<double> xs{0,0,0,0};
vector<double> selectX;
vector<double> selectY;
selectX.push_back(..);
...
selectY.push_back(..);
...
double R = polynomialFit(selectX,selectY,3,xs);