如何在C++ Eigen中复现Python statsmodels.api.OLS的回归结果？

我用Python的statsmodels计算两个向量的回归直线斜率和截距，代码如下：

A = [1,2,5,7,14,17,19]
b = [2,14,6,7,13,27,29]

A = sm.add_constant(A)
results = sm.OLS(A, b).fit()

print("results: ", results.params)

输出: [0.04841897 0.64278656]

现在想用C的Eigen库复现这个结果，但添加全1列后的回归结果和Python不一致，不加列或加全0列时斜率和Python一致但截距为0。我的C代码：

Eigen::VectorXd A(7);
Eigen::VectorXd b(7);

A << 1,2,5,7,14,17,19;
b << 2,14,6,7,13,27,29;

MatrixXd new_A(A.rows(), 2);
VectorXd d = VectorXd::Constant(A.rows(), 1);
new_A << A, d;

Eigen::MatrixXd res = new_A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);

cout << " slope: " << res.coeff(0, 0) << " intercept: " << res.coeff(1, 0) << endl;
cout << "dbl check: " << (new_A.transpose() * new_A).ldlt().solve(new_A.transpose() * b) << endl;

添加全1列的输出 -> slope: 1.21644 intercept:2.70444；不加列或加全0列的输出 -> slope: 0.642787 intercept: 0

想问：如何在C++中得到和Python一致的结果？哪组结果是正确的？我更倾向信任Python的结果，因为不加列时斜率和Python一致。

问题根源：OLS参数顺序搞反了 #

Python代码里有个容易踩的坑：sm.OLS(A, b).fit()的参数顺序写反了。statsmodels的OLS类第一个参数是因变量(y)，第二个是自变量(X)，你把加了常数项的A放在第一个位置，b放第二个，相当于拟合的是A = b*β₀ + 1*β₁，而非常规的b = β₀*1 + β₁*A。这就解释了为什么不加列时两边斜率一致——此时都是拟合无截距的b = β*A，参数顺序不影响结果，但加截距后逻辑差异就显现了。

正确处理方式 #

1. 修正Python代码（得到真实回归结果） #

如果你的需求是拟合常规一元线性回归y = β₀ + β₁x（x为A，y为b），Python代码应该改成：

import statsmodels.api as sm

A = [1,2,5,7,14,17,19]
b = [2,14,6,7,13,27,29]

# 自变量矩阵X：添加常数项列对应截距β₀
X = sm.add_constant(A)
# OLS(因变量, 自变量)
results = sm.OLS(b, X).fit()

print("results: ", results.params)

输出: [2.70443686 1.21643836]
这个结果和你C中添加全1列的输出完全一致，说明C的代码逻辑是对的，原Python代码只是参数顺序错了。

2. 复现原Python错误结果的C++代码 #

如果你非要复现原Python的反向拟合结果（即A = β₀*b + β₁*1），只需要调换Eigen中的求解对象和自变量矩阵：

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

using namespace Eigen;
using namespace std;

int main() {
    Eigen::VectorXd A(7);
    Eigen::VectorXd b(7);

    A << 1,2,5,7,14,17,19;
    b << 2,14,6,7,13,27,29;

    // 构造自变量矩阵：b作为x，添加全1列
    MatrixXd new_b(b.rows(), 2);
    VectorXd d = VectorXd::Constant(b.rows(), 1);
    new_b << b, d;

    // 拟合A = β₀*b + β₁*1
    Eigen::VectorXd res = new_b.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(A);

    cout << " 对应Python原结果的参数: [" << res(1) << ", " << res(0) << "]" << endl;
    // 输出会接近 [0.04841897, 0.64278656]
    return 0;
}

注意：原Python输出的[0.0484, 0.6428]对应A = 0.6428*b + 0.0484*1，所以C++中res(0)是b的系数，res(1)是常数项，和Python输出顺序相反。

结论 #

• 正确的常规回归结果是C++添加全1列的输出：截距2.7044，斜率1.2164，对应拟合式b = 2.7044 + 1.2164*A。
• 原Python代码因参数顺序错误，得到的是反向拟合的结果，并非你需要的回归直线。

内容的提问来源于stack exchange，提问作者Merlin