【基础学习】拟合点集及拟合求导

yoyo-han

已于 2024-04-28 13:30:52 修改

阅读量167

点赞数 1

分类专栏：无人机文章标签：最小二乘法算法学习

于 2024-04-15 12:45:44 首次发布

本文链接：https://blog.csdn.net/weixin_45058446/article/details/137776342

版权

无人机专栏收录该内容

4 篇文章 0 订阅

订阅专栏

主要是为了手撕题准备

直线拟合

直线拟合是假设模型为一条直线 y=kx+b，使得观测数据点到该直线的距离的平方和fx最小
原理如图：
在这里插入图片描述

#include <iostream>
#include <vector>

using namespace std;
                                                                                                                         
void leastSquaresFit(const vector<double>& x, const vector<double>& y, double& k, double& b) {
    int n = x.size();
    double sum_x = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_x2 = 0.0;

    for (int i = 0; i < n; i++) {
        sum_x += x[i];
        sum_y += y[i];
        sum_xy += x[i] * y[i]; 
        sum_x2 += x[i] * x[i];
    }

    k = (sum_xy * n - sum_x * sum_y) / (sum_x2 * n - sum_x * sum_x);
    b = (sum_y * sum_x2 - sum_x * sum_xy) / (sum_x2 * n - sum_x * sum_x);
}

int main() {
    vector<double> x = {1, 2, 3, 4, 5};
    vector<double> y = {2, 3, 4, 5, 6};
    double k, b;
    leastSquaresFit(x, y, k, b);
    cout << "y = " << k << "x + " << b << endl;
    return 0;
}

线性回归

线性回归对应 fx = a0+a1x1+a2x2……an-1xn-1的情况，其中x相当于特征量，要使||fx - y||尽可能小，原理如图
在这里插入图片描述
不调用库的形式

#include <iostream>
#include <vector>

using namespace std;

void linearRegression(vector<vector<double>>& X, vector<double>& y, vector<double>& coefficients) {
    int n = X.size(); 
    int m = X[0].size(); 
    vector<vector<double>> X_transpose(m, vector<double>(n)); 

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            X_transpose[j][i] = X[i][j];
        }
    }

    vector<vector<double>> X_transpose_X(m, vector<double>(m));
    for (int i = 0; i < m;  i++) {
        for (int j = 0; j < m; j++) {
            double sum = 0.0;
            for (int k = 0; k < n; k++) {
                sum += X_transpose[i][k] * X[k][j];
            }
            X_transpose_X[i][j] = sum;
        }
    }

    vector<double> X_transpose_y(m);
    for (int i = 0; i < m; i++) {
        double sum = 0.0;
        for (int j = 0; j < n; j++) {
            sum += X_transpose[i][j] * y[j];
        }
        X_transpose_y[i] = sum;
    }

    // 调整矩阵为上三角形式 求解coefficients(列向量）
    for (int i = 0; i < m; i++) {
        double divisor = X_transpose_X[i][i];
        for (int j = i; j < m; j++) {
            X_transpose_X[i][j] /= divisor; // 对角线变为1 仅对上三角操作
        }
        X_transpose_y[i] /= divisor;
        for (int k = i + 1; k < m; k++) {
            double factor = X_transpose_X[k][i];
            for (int j = i; j < m; j++) {
                X_transpose_X[k][j] -= factor * X_transpose_X[i][j]; // 逐列消掉下三角
            }
            X_transpose_y[k] -= factor * X_transpose_y[i];
        }
    }
    // 得到上三角形式
    for (int i = m - 1; i >= 0; i--) {
        coefficients[i] = X_transpose_y[i];
        for (int j = i + 1; j < m; j++) {
            coefficients[i] -= X_transpose_X[i][j] * coefficients[j];
        }
    }
}

int main() {
    // 默认存在x0=1
    vector<vector<double>> X = {{1, 1,4}, {1, 2,5}, {1, 3,3}, {1, 4,2}, {1, 5,1}};
    vector<double> y = {6,8,7,7,7};
    vector<double> coefficients(X[0].size()); 
    linearRegression(X, y, coefficients);

    cout << "y = " << coefficients[0];
    for (int i = 1; i < coefficients.size(); i++) {
        cout << " +";
        cout << coefficients[i] << "x" << i;
    }
    return 0;
}

使用Eigen库的情况

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;

void linearRegression(const MatrixXd& X, const VectorXd& y, VectorXd& coefficients) {
    coefficients = (X.transpose() * X).inverse() * X.transpose() * y;
}

int main() {
    MatrixXd X(5, 2); 
    VectorXd y(5); 
    X << 1, 1,
         1, 2,
         1, 3,
         1, 4,
         1, 5,
    y << 2, 3, 4, 5, 6;
    VectorXd coefficients(2); 
    linearRegression(X, y, coefficients);
    cout << "y = " << coefficients[0];
    for (int i = 1; i < coefficients.size(); ++i) {
        cout << " +";
        cout << coefficients[i] << "x" << i;
    }
    return 0;
}

多项式拟合

形如fx =a0+a1x¹+a2x²+……的形式
用库函数版本与线性拟合基本类似，只是x矩阵的格式不同，纵轴为degree，但求coefficients的公式仍然适用

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

void polynomialFit(vector<double>& x, vector<double>& y, int degree, vector<double>& coefficients) {
    int n = x.size();
    int m = degree + 1;
    vector<vector<double>> X(n, vector<double>(m));

    // 构造自变量矩阵
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            X[i][j] = pow(x[i], j);
        }
    }

    // 构造增广矩阵，增广矩阵的最后一列是因变量向量
    vector<vector<double>> augmentedX(n, vector<double>(m + 1));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            augmentedX[i][j] = X[i][j];
        }
        augmentedX[i][m] = y[i];
    }


    // 使用高斯消元法求解增广矩阵
    for (int i = 0; i < m; ++i) {
        for (int j = i + 1; j < n; ++j) {
            double ratio = augmentedX[j][i] / augmentedX[i][i];
            for (int k = i; k < m + 1; ++k) {
                augmentedX[j][k] -= ratio * augmentedX[i][k];
            }
        }
    }

    // 回代求解未知系数
    coefficients.resize(m);
    for (int i = m - 1; i >= 0; --i) {
        coefficients[i] = augmentedX[i][m];
        for (int j = i + 1; j < m; ++j) {
            coefficients[i] -= augmentedX[i][j] * coefficients[j];
        }
        coefficients[i] /= augmentedX[i][i];
    }

}

int main() {
    vector<double> x = {1, 2, 3, 4, 5}; // 自变量向量
    vector<double> y = {1, 4, 9, 16, 25}; // 因变量向量
    int degree = 3; // 多项式的阶数
    vector<double> coefficients; 
    polynomialFit(x, y, degree, coefficients);
    cout << "Polynomial Fit Model: fx = " << coefficients[0];
    for (int i = 1; i < coefficients.size(); ++i) {
        cout << " + ";
        cout << coefficients[i] << "x^" << i;
    }
    cout << endl;
    return 0;
}

多项式求导

输入多项式的最大阶次和求导次数，自变量的值及系数，输出求导后的多项式在自变量 x 处的值

// 输入多项式的最大阶次和求导次数，自变量的值及系数，输出求导后的多项式在 x 处的值
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

void caculate_derivative(vector<int>& coeff, int m, int n){
    for(int i=0;i<m;i++){
        for (int j = 0; j <= n; j++)
        {
            coeff[j] *= (j-i) ? (j-i) : 0;
        }
        
    }
}

int evaluate(vector<int> &coeff, int x, int m){
    int res=0;
    for(int i=0;i<coeff.size();++i){
        if(i<m) continue;
        else if(i==m) res = coeff[i];
        else{
            res += coeff[i] * pow(x, (i-m));
        }
    }
    return res;
}

int main(){
    int n, m, x;
    cin >> n >> m >> x;
    vector<int> coeff(n+1);
    for(int i=0;i<=n;i++){
        cin >> coeff[i];
    }
    caculate_derivative(coeff, m, n);
    int res = evaluate(coeff, x, m);
    cout << res << endl;
    return 0;
}