使用 Armadillo 在 C++ 中实现 numpy.polyfit 和 numpy.polyval

Implementation of numpy.polyfit and numpy.polyval in C++ with Armadillo

我正在尝试重现 numpy.polyfit and the following application of numpy.polyval in C++ by using Armadillo 的结果。

这是我的尝试:

using namespace arma;

vec fastLm(const vec& y,
           const mat& X,
           int order)
{
  mat extended_X(X);
  
  // Column bind the higher order regressors to the initial matrix
  for(int i = 2; i < order + 1; ++i)
  {
    extended_X = join_rows(extended_X, pow(X, i));
  }
  
  // Join another column, made of '1', in order to have the intercept
  extended_X = join_rows(mat(X.n_rows, 1, fill::ones), extended_X);
  
  // Solve the linear regression by OLS
  colvec coef = solve(extended_X, y);

  // Return the fit
  return extended_X * coef;
}

我希望得到与以下相同的结果:

import numpy as np

def fastLm(y, X, order):

    # Fit the polynomial regression
    rg = np.polyfit(X, y, order)

    # Return the fit
    C = np.polyval(rg, X)

    return C

但是,我的测试显示出差异和奇怪的结果,我很难找到和调试其原因。你能告诉我我的“翻译”是否正确或修复它吗?

对我来说一切都很好,已经用

试过你的功能
int main()
{
  mat x=linspace(0,1,5);
  vec y=1/(1+x);
  y.print("Y");
  mat Yhat = fastLm(y,x,3);
  Yhat.print("Yhat");
}

给出结果

Y
   1.0000
   0.8000
   0.6667
   0.5714
   0.5000
Yhat
   0.9998
   0.8008
   0.6654
   0.5722
   0.4998

你的python代码对应的结果是

[1.         0.8        0.66666667 0.57142857 0.5       ]
[0.99979592 0.80081633 0.66544218 0.5722449  0.49979592]

... 和 Matlab

>> Y
Y =

   1.00000   0.80000   0.66667   0.57143   0.50000

>> Yhat
Yhat =

   0.99980   0.80082   0.66544   0.57224   0.49980

>>