带截距的多元分位数回归 C++带类完整的正确实现及成功案例

下面是一个带截距的多元分位数回归的C++类的完整实现及一个成功的案例：

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

#include "lp_solve/lp_lib.h"

class QuantileRegression {
public:
    QuantileRegression(double q) : tau(q) {}

    // 添加一个样本，x是自变量向量，y是因变量
    void addSample(const std::vector<double>& x, double y) {
        samples.emplace_back(x, y);
    }

    // 训练模型，返回估计的系数向量
    std::vector<double> train() {
        int n = samples.size();
        int m = samples[0].first.size();

        // 构造优化问题的矩阵形式
        std::vector<std::vector<double>> A(n + 2 * m, std::vector<double>(n + m, 0.0));
        std::vector<double> b(n + 2 * m, 0.0);
        for (int i = 0; i < n; i++) {
            b[i] = samples[i].second;
            for (int j = 0; j < m; j++) {
                A[i][j] = samples[i].first[j];
            }
        }

        // 添加截距项
        for (int i = 0; i < n; i++) {
            A[i][m + i] = 1.0;
        }

        // 添加约束
        for (int i = 0; i < m; i++) {
            A[n + i][i] = 1.0;
            A[n + m + i][n + i] = 1.0;
            b[n + i] = tau * n;
            b[n + m + i] = (tau - 1) * n;
        }

        // 求解线性规划问题
        std::vector<double> c(n + m, 0.0);
        std::fill(c.begin(), c.begin() + n, 0.0);
        std::fill(c.begin() + n, c.end(), 1.0);

        std::vector<double> x(n + m, 0.0);
        lp_solve(n + m, A, b, c, x);

        // 提取系数向量
        std::vector<double> beta(m + 1, 0.0);
        for (int i = 0; i <= m; i++) {
            beta[i] = x[i];
        }
        return beta;
    }

private:
    double tau;
    std::vector<std::pair<std::vector<double>, double>> samples;

    // 使用线性规划库求解问题
    void lp_solve(int n, const std::vector<std::vector<double>>& A,
                  const std::vector<double>& b, const std::vector<double>& c,
                  std::vector<double>& x) {
        lprec* lp = make_lp(0, n);
        set_verbose(lp, IMPORTANT);
        set_add_rowmode(lp, TRUE);

        std::vector<int> colno(n + 1);
        for (int i = 0; i <= n; i++) {
            colno[i] = i + 1;
        }

        for (int i = 0; i < A.size(); i++) {
            add_constraintex(lp, A[i].size(), &A[i][0], &colno[0], GE, b[i]);
        }

        set_obj_fnex(lp, c.size(), &c[0], &colno[0]);

        set_add_rowmode(lp, FALSE);
        set_minim(lp);

        int ret = solve(lp);
        if (ret == OPTIMAL) {
            get_variables(lp, &x[0]);
        } else {
            std::cerr << "LP solve failed: " << ret << std::endl;
        }

        delete_lp(lp);
    }
};

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 = {1.0, 2.0, 3.0, 4.0};

    // 训练模型
    QuantileRegression model(0.5);
    for (int i = 0; i < X.size(); i++) {
        model.addSample(X[i], y[i]);
    }
    std::vector<double> beta = model.train();

    // 打印结果
    std::cout << "Coefficients: ";
    for (int i = 0; i < beta.size(); i++) {
        std::cout << beta[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}

这里使用了一个开源的线性规划库 lp_solve 来求解问题。在运行程序前需要先安装该库并将头文件和库文件正确地链接到程序中。在这个例子中，我们使用了一个包含四个样本的数据集，其中每个样本有两个自变量和一个因变量。我们使用了分位数 q=0.5 计算的结果，即中位数回归。运行程序后，我们将获得估计的系数向量，其包含两个元素，分别对应于两个自变量和一个截距项。