Реализация numpy .polyfit и numpy .polyval на C ++ с помощью Armadillo - PullRequest
Новая
       276

Реализация numpy .polyfit и numpy .polyval на C ++ с помощью Armadillo

0 голосов
/ 07 августа 2020

Я пытаюсь воспроизвести результаты numpy.polyfit и следующего приложения numpy.polyval на C ++, используя 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

Однако мои тесты показывают несоответствия и странные результаты, причина которых я m возникают трудности с поиском и отладкой. Подскажите, пожалуйста, верен ли мой "перевод" или поправьте?

1 Ответ

1 голос
/ 08 августа 2020

Для меня все выглядит хорошо, пробовал вашу функцию с 

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

>>
...