回归模型——稳健回归 C++带类实现及案例

稳健回归（Robust Regression）是一种针对数据异常值影响的回归分析方法。相比于传统的最小二乘法，稳健回归能够忽略掉数据中的异常值，从而得到更加可靠的回归结果。

C++中可以使用类来实现稳健回归算法，主要包括以下步骤：

1. 定义类RobustRegression，包含以下成员变量和成员函数：

成员变量：
- vector<double> x：自变量
- vector<double> y：因变量
- double b0：截距
- double b1：斜率
- int n：样本量
- double w：迭代权重

成员函数：
- 构造函数RobustRegression(vector<double>& x, vector<double>& y)：初始化自变量和因变量
- double getB0()：获取截距
- double getB1()：获取斜率
- double getR()：获取相关系数
- void fit()：拟合回归模型

2. 在构造函数中初始化自变量和因变量。在拟合回归模型时，使用Huber估计量来计算迭代权重。具体实现如下：

void RobustRegression::fit() {
    double delta = 1.345 * median(abs(x - median(x))); // Huber阈值
    vector<double> residuals(n); // 残差
    double w_sum = 0; // 权重和
    for (int i = 0; i < n; i++) {
        residuals[i] = y[i] - (b0 + b1 * x[i]); // 计算残差
        if (abs(residuals[i]) <= delta) {
            w = 1; // 权重为1
        } else {
            w = delta / abs(residuals[i]); // 权重为delta/|residuals|
        }
        w_sum += w; // 累加权重
        b0 += w * residuals[i]; // 更新截距
        b1 += w * residuals[i] * x[i]; // 更新斜率
    }
    b0 /= w_sum; // 计算截距的加权平均值
    b1 /= w_sum; // 计算斜率的加权平均值
}

3. 在类中实现获取截距、斜率和相关系数的成员函数。其中相关系数可以使用Pearson相关系数来计算：

double RobustRegression::getR() {
    double x_sum = 0, y_sum = 0, xy_sum = 0, x2_sum = 0, y2_sum = 0;
    for (int i = 0; i < n; i++) {
        x_sum += x[i];
        y_sum += y[i];
        xy_sum += x[i] * y[i];
        x2_sum += x[i] * x[i];
        y2_sum += y[i] * y[i];
    }
    double r = (n * xy_sum - x_sum * y_sum) / sqrt((n * x2_sum - x_sum * x_sum) * (n * y2_sum - y_sum * y_sum));
    return r;
}

4. 编写一个测试程序，读入自变量和因变量数据，创建RobustRegression对象并拟合回归模型，输出截距、斜率和相关系数。

完整的C++类实现代码如下：

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

using namespace std;

class RobustRegression {
public:
    RobustRegression(vector<double>& x, vector<double>& y) : x(x), y(y), b0(0), b1(0), n(x.size()), w(0) {}
    double getB0() { return b0; }
    double getB1() { return b1; }
    double getR();
    void fit();
private:
    vector<double> x;
    vector<double> y;
    double b0;
    double b1;
    int n;
    double w;
};

void RobustRegression::fit() {
    double delta = 1.345 * median(abs(x - median(x))); // Huber阈值
    vector<double> residuals(n); // 残差
    double w_sum = 0; // 权重和
    for (int i = 0; i < n; i++) {
        residuals[i] = y[i] - (b0 + b1 * x[i]); // 计算残差
        if (abs(residuals[i]) <= delta) {
            w = 1; // 权重为1
        } else {
            w = delta / abs(residuals[i]); // 权重为delta/|residuals|
        }
        w_sum += w; // 累加权重
        b0 += w * residuals[i]; // 更新截距
        b1 += w * residuals[i] * x[i]; // 更新斜率
    }
    b0 /= w_sum; // 计算截距的加权平均值
    b1 /= w_sum; // 计算斜率的加权平均值
}

double RobustRegression::getR() {
    double x_sum = 0, y_sum = 0, xy_sum = 0, x2_sum = 0, y2_sum = 0;
    for (int i = 0; i < n; i++) {
        x_sum += x[i];
        y_sum += y[i];
        xy_sum += x[i] * y[i];
        x2_sum += x[i] * x[i];
        y2_sum += y[i] * y[i];
    }
    double r = (n * xy_sum - x_sum * y_sum) / sqrt((n * x2_sum - x_sum * x_sum) * (n * y2_sum - y_sum * y_sum));
    return r;
}

int main() {
    vector<double> x = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    vector<double> y = {2.1, 4.2, 6.3, 8.4, 10.5, 12.6, 14.7, 16.8, 18.9, 21};
    RobustRegression rr(x, y);
    rr.fit();
    cout << "b0 = " << rr.getB0() << endl;
    cout << "b1 = " << rr.getB1() << endl;
    cout << "r = " << rr.getR() << endl;
    return 0;
}

上述程序中，自变量和因变量数据都是手动输入的。在实际应用中，可以从文件或数据库中读取数据。