基于c++与eigen库的最小二乘法

最新推荐文章于 2024-04-03 14:08:36 发布

苏州安瑞软件有限公司

最新推荐文章于 2024-04-03 14:08:36 发布

阅读量247

点赞数 1

分类专栏：数学计算文章标签： c++ 最小二乘法

本文链接：https://blog.csdn.net/qq_37033647/article/details/131248778

版权

数学计算专栏收录该内容

3 篇文章 1 订阅

订阅专栏

该代码实现了一个C++类FitCurve，利用Eigen库进行二次多项式曲线拟合。它接受散点坐标值列表，计算多项式系数并返回拟合的方差R²。Vandermonde矩阵被用于求解系数，并通过比较残差平方和（SSR）与总平方和（SST）来计算R²。

摘要由CSDN通过智能技术生成

一、涉及的工具和库

我编译用的工具是Qt，不过本文中未设计任何Qt相关的类，用的是C++和eigen库。

二、头文件

#ifndef FITCURVE_H
#define FITCURVE_H

#include "FitCurve_global.h"
#include "vector"
#include "numeric"
#include "Eigen/Dense"
using namespace std;
class FitCurve
{
public:
    FitCurve();
    double quadraticFitCurve(vector<float> &X,vector<float> 
};
#endif // FITCURVE_H

三、源文件

函数的X，Y为散点坐标值列表，orders代表多项式最高次数，coef为存放多项式系数的数组，coef数组长度为orders+1
返回值是方差R²

#include "fitcurve.h"

FitCurve::FitCurve()
{
}
double quadraticFitCurve(vector<float> &X, vector<float> &Y, uint8_t orders, double *coef)
{
    if (X.size() < 2 || Y.size() < 2 || X.size() != Y.size() || orders < 1)
            exit(EXIT_FAILURE);
    // map sample data from STL vector to eigen vector
    Eigen::Map<Eigen::VectorXf> sampleX(X.data(), X.size());
    Eigen::Map<Eigen::VectorXf> sampleY(Y.data(), Y.size());
    Eigen::MatrixXf mtxVandermonde(X.size(), orders + 1);  // Vandermonde matrix of X-axis coordinate vector of sample data
    Eigen::VectorXf colVandermonde = sampleX;              // Vandermonde column

    // construct Vandermonde matrix column by column
    for (size_t i = 0; i < orders + 1; ++i)
    {
        if (0 == i)
        {
            mtxVandermonde.col(0) = Eigen::VectorXf::Constant(X.size(), 1, 1);
            continue;
        }
        if (1 == i)
        {
            mtxVandermonde.col(1) = colVandermonde;
            continue;
        }
        colVandermonde = colVandermonde.array()*sampleX.array();
        mtxVandermonde.col(i) = colVandermonde;
    }
    // calculate coefficients vector of fitted polynomial
    Eigen::VectorXf result = (mtxVandermonde.transpose()*mtxVandermonde).inverse()*(mtxVandermonde.transpose())*sampleY;
    for(uint8_t j = 0; j < orders + 1;j++){
        coef[j]=result[j];
    }

    float sum = std::accumulate(std::begin(Y),std::end(Y),0.0);
    float mean = sum/Y.size();//求Y的平均值
    double SSR=0;
    double SST=0;
    for(uint8_t j=0;j<Y.size();j++){
        float yi = 0;
        for(uint8_t k=0;k<orders+1;k++){
            yi += result[k]*pow(X[j],k);
        }
        SSR += pow(yi-Y[j],2);
        SST += pow(Y[j]-mean,2);
    }
    double R = 1-SSR/SST;
    return R;
}

四、调用

#include <vector>
vector<float> x,y;
x.push_back(192);
x.push_back(190);
x.push_back(187);
x.push_back(184);
x.push_back(178);
x.push_back(170);
x.push_back(163);
x.push_back(153);
x.push_back(138);
x.push_back(120);
y.push_back(191);
y.push_back(189);
y.push_back(187);
y.push_back(185);
y.push_back(181);
y.push_back(175);
y.push_back(170);
y.push_back(163);
y.push_back(153);
y.push_back(139);
double xs[3];
double R=quadraticFitCurve(x,y,2,xs);