带截距的多元稳健回归 C++ 实现及案例

多元稳健回归（Robust Regression）是一种针对数据中存在异常值或者噪声的回归分析方法。它的主要特点在于对异常值的抗干扰能力比较强。下面是带截距的多元稳健回归的 C++ 实现及一个简单的案例。

首先，需要实现一个计算加权样本均值和方差的函数：

#include <vector>
#include <cmath>

void weighted_mean_var(const std::vector<double>& x, const std::vector<double>& w, double& mean, double& var) {
    double sum_w = 0.0;
    double sum_wx = 0.0;
    double sum_wx2 = 0.0;
    for (int i = 0; i < x.size(); i++) {
        sum_w += w[i];
        sum_wx += w[i] * x[i];
        sum_wx2 += w[i] * x[i] * x[i];
    }
    mean = sum_wx / sum_w;
    var = (sum_wx2 / sum_w) - std::pow(mean, 2.0);
}

然后，实现带截距的多元稳健回归模型：

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

class RobustRegression {
public:
    RobustRegression(double tol=1e-4, int max_iter=1000) : tol_(tol), max_iter_(max_iter) {}
    void fit(const std::vector<std::vector<double>>& X, const std::vector<double>& y);
    double predict(const std::vector<double>& x);

private:
    double tol_;
    int max_iter_;
    Eigen::VectorXd coef_;
    Eigen::VectorXd resid_;
    Eigen::VectorXd weights_;

    void update_weights();
    void update_coef();
};

void RobustRegression::fit(const std::vector<std::vector<double>>& X, const std::vector<double>& y) {
    int n = X.size();
    int p = X[0].size();

    coef_ = Eigen::VectorXd::Zero(p + 1);
    resid_ = Eigen::VectorXd::Zero(n);
    weights_ = Eigen::VectorXd::Ones(n);

    Eigen::MatrixXd X_mat(n, p + 1);
    for (int i = 0; i < n; i++) {
        X_mat(i, 0) = 1.0;
        for (int j = 0; j < p; j++) {
            X_mat(i, j + 1) = X[i][j];
        }
    }

    int iter = 0;
    while (iter < max_iter_) {
        Eigen::VectorXd y_pred = X_mat * coef_;
        resid_ = y - y_pred;
        update_weights();
        update_coef();
        iter++;
    }
}

double RobustRegression::predict(const std::vector<double>& x) {
    Eigen::VectorXd x_vec(x.size() + 1);
    x_vec << 1.0, Eigen::VectorXd::Map(&x[0], x.size());
    return x_vec.dot(coef_);
}

void RobustRegression::update_weights() {
    double median_resid = resid_.median();
    double mad_resid = (resid_.array() - median_resid).array().abs().median();
    double c = 1.4826 * mad_resid;
    weights_ = (resid_.array() / c).array().square().array().max(1.0).matrix();
}

void RobustRegression::update_coef() {
    Eigen::MatrixXd W_mat = weights_.asDiagonal();
    Eigen::MatrixXd XtW = X_mat_.transpose() * W_mat;
    Eigen::MatrixXd XtWX = XtW * X_mat_;
    Eigen::VectorXd XtWy = XtW * y_;
    coef_ = XtWX.colPivHouseholderQr().solve(XtWy);
}

最后，给出一个简单的案例：

#include <iostream>
#include "robust_regression.h"

int main() {
    std::vector<std::vector<double>> X = {{1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0}, {4.0, 5.0}};
    std::vector<double> y = {2.0, 3.0, 4.0, 5.0};

    RobustRegression model;
    model.fit(X, y);

    std::vector<double> x = {5.0, 6.0};
    double y_pred = model.predict(x);
    std::cout << "Predicted value: " << y_pred << std::endl;

    return 0;
}