用犰狳在 C++ 中实现 numpy.polyfit 和 numpy.polyval

Posted 2023-02-17

【中文标题】用犰狳在 C++ 中实现 numpy.polyfit 和 numpy.polyval

【英文标题】：Implementation of numpy.polyfit and numpy.polyval in C++ with Armadillo

【发布时间】：2020-08-07 09:16:16

【问题描述】：

我正在尝试使用 Armadillo 在 C++ 中重现 numpy.polyfit 的结果以及 numpy.polyval 的以下应用。

这是我的尝试：

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

但是，我的测试显示出差异和奇怪的结果，导致我难以查找和调试。您能否告诉我我的“翻译”是否正确或修复它？

【问题讨论】：

为什么要重新实现？犰狳已经有了polyfit 和polyval 函数。

真的吗？不知道:(

【参考方案1】：

对我来说一切都很好，已经尝试过你的功能

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

>>