用C++实现迭代重加权最小二乘法（IRLS）

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

using namespace Eigen;

VectorXd irls(const VectorXd& x, const VectorXd& y, int degree, int iterations, double epsilon) {
    int m = x.size();
    MatrixXd X(m, degree + 1);
    
    // 构建特征矩阵X
    for (int i = 0; i < m; i++) {
        for (int j = 0; j <= degree; j++) {
            X(i, j) = pow(x(i), j);
        }
    }
    
    VectorXd theta = VectorXd::Zero(degree + 1);
    
    // 迭代更新参数
    for (int iter = 0; iter < iterations; iter++) {
        VectorXd y_pred = X * theta;
        VectorXd residuals = y - y_pred;
        MatrixXd weights = MatrixXd::Zero(m, m);
        
        // 计算权重矩阵
        for (int i = 0; i < m; i++) {
            weights(i, i) = 1 / (fabs(residuals(i)) + epsilon);
        }
        
        // 更新参数
        theta = (X.transpose() * weights * X).inverse() * X.transpose() * weights * y;
    }
    
    return theta;
}

int main() {
    // 生成样本数据
    int n = 100;
    VectorXd x(n);
    VectorXd y(n);
    for (int i = 0; i < n; i++) {
        x(i) = i / 10.0;
        y(i) = 2 * x(i) + 1 + (rand() / (RAND_MAX + 1.0)) - 0.5;
    }
    
    // 使用IRLS拟合曲线
    int degree = 1;  // 多项式的次数
    int iterations = 10;  // 迭代次数
    double epsilon = 0.01;  // 避免除零错误的小值
    VectorXd theta = irls(x, y, degree, iterations, epsilon);
    
    // 输出拟合参数
    std::cout << "拟合参数：" << std::endl;
    for (int i = 0; i <= degree; i++) {
        std::cout << "theta[" << i << "] = " << theta(i) << std::endl;
    }
    
    return 0;
}

在上述代码中，我们使用了Eigen库来进行矩阵和向量的计算。首先，我们定义了一个irls函数，该函数接受输入的x和y向量，以及多项式的次数、迭代次数和epsilon（用于避免除零错误的小值）作为参数。函数内部首先构建特征矩阵X，然后初始化参数theta为全零向量。接下来，使用迭代的方式更新参数theta，直到达到指定的迭代次数。在每次迭代中，计算预测值y_pred、残差residuals和权重矩阵weights，并使用加权最小二乘法更新参数theta。最后，返回更新后的参数theta。