c++ - c++中具有奇数和偶数指数的多项式函数

Question

我是新来的，也是编程的初学者。
对于一个小型的私有(private)使用项目，我想计算各种多项式模型的系数。一种数据集的 polyfit。

正态多项式模型 (x^0, x^1, x^2, x^3)

只有偶数指数 (x^0, x^2, x^4)

只有奇数指数

我在下一页找到了一个应用正态多项式模型的函数(这段代码不是我写的，权利不是我的)
https://gist.github.com/chrisengelsma/108f7ab0a746323beaaf7d6634cf4add
但我想对此代码进行建模以使用其他两个模型。
不幸的是，我的数学技能有限。
如果我看到正确的，我只需要更改以下几行或？
原来的:

     // X = vector that stores values of sigma(xi^2n)
  std::vector<TYPE> X(tnp1);
  for (int i = 0; i < tnp1; ++i) {
    X[i] = 0;
    for (int j = 0; j < N; ++j)
      X[i] += (TYPE)pow(x.get(j), i);
  }

正态多项式的原始值。
下一个偶数

     // X = vector that stores values of sigma(xi^2n)
  std::vector<TYPE> X(tnp1);
  for (int i = 0; i < tnp1; i+2) {
    X[i] = 0;
    for (int j = 0; j < N; ++j)
      X[i] += (TYPE)pow(x.get(j), i);
  }

下一个是奇数

     // X = vector that stores values of sigma(xi^2n)
  std::vector<TYPE> X(tnp1);
  for (int i = 1; i < tnp1; i+2) {
    X[i] = 0;
    for (int j = 0; j < N; ++j)
      X[i] += (TYPE)pow(x.get(j), i);
  }

这个片段也是一样的:

  // Y = vector to store values of sigma(xi^n * yi)
  std::vector<TYPE> Y(np1);
  for (int i = 0; i < np1; ++i) {
    Y[i] = (TYPE)0;
    for (int j = 0; j < N; ++j) {
      Y[i] += (TYPE)pow(x[j], i)*y[j];
    }
  }

我无法验证结果。至少我不知道怎么做。
因此，我希望有人能告诉我，我已经正确理解了这一点并可以实现。

Answer 1

如果您的拟合函数相对于估计系数是线性的，则可以使用线性最小二乘算法作为一般方法。我将使用 alglib对于示例，但任何 C++ 数学库也应该完成这项工作。
您的拟合函数可以表示为

f(t) = c0 * f0(t) + c1 * f1(t) + c2 * f2(t) + ...

在哪里(在您的偶数多项式示例中)您的 fx 可以是

f0 = constant
f1 = t^2
f2 = t^4

有了它，您可以定义 Least Squares问题:

A * c - y  ---> min

在哪里 A是您在不同时间点拟合 f(t) 曲线值的矩阵

A = [[f0(t0), f1(t0), f2(t0)],
  [f0(t1), f1(t1), f2(t1)]
  [f0(t2), f1(t2), f2(t2)],
          ...
  [f0(tn), f1(tn), f2(tn)]];

c是要估计的系数 vector ， y是您测量的值( y 和 t 是您对算法的输入)。
这个问题可以通过 alglib::lsfitlinear 来解决:

    const size_t numPoints(1000);
    double p0(-1000);
    double p2(-5.5);
    double p4(1.1);
    std::vector<double> t(numPoints);
    std::generate(t.begin(), t.end(), [i = 0]() mutable { return 0.001*i++; });
    std::random_device rd{};
    std::mt19937 gen{rd()};
    std::normal_distribution<> gauss{0, 0.5};  // let's apply some noise
    std::vector<double> y(numPoints);
    // generate fake data
    std::generate(y.begin(), y.end(), [i = 0, &t, &gauss, &gen, p0, p2, p4]() mutable {
        double result = p0 + p2*std::pow(t[i], 2) + p4*std::pow(t[i], 4);
        result += gauss(gen);
        i++;
        return result;
    });
    alglib::real_1d_array y_array;
    alglib::real_1d_array t_array;
    alglib::real_2d_array F_array;
    y_array.setcontent(y.size(), y.data());
    t_array.setcontent(t.size(), t.data());
    // LS array for 3 basis functions: f0 = 1, f1 = x^2, f2 = x^4
    F_array.setlength(numPoints, 3);
    std::normal_distribution<> gauss2{1, 1}; 
    for (size_t i = 0; i < numPoints; i++)
    {
        F_array(i, 0) = 1;
        F_array(i, 1) = std::pow(t[i], 2);
        F_array(i, 2) = std::pow(t[i], 4);
    }
    alglib::ae_int_t info;
    alglib::real_1d_array coeffs;
    alglib::lsfitreport rep;
    alglib::lsfitlinear(y_array, F_array, numPoints, 3, info, coeffs, rep);
    std::cout << "Success: " << info << std::endl;
    std::cout << "Results: " << coeffs(0) << ", " << coeffs(1) << ", " << coeffs(2) << std::endl;
    std::cout << "RMS: " << rep.rmserror << std::endl;

但是，对于较大的值，这种方法可能会在数值上变得不稳定。尝试设置更大的时间值，看看结果会比下面的替代方法更差。请注意，您的矩阵是 Vandermonde matrix ，因此时间值分布在 Y 轴两侧将变得不可逆，因为您只使用偶数指数系数。示例结果:

Success: 1
Results: -1000.04, -5.45533, 1.19111
RMS: 0.479453

幸运的是，在这种特殊情况下(甚至多项式系数)，您可以通过替换来完成数学技巧

tp <--- t^2

并将您的问题改造成具有函数的标准多项式拟合

f(t) = C + t^2 + t^4 = c + tp + tp^2

并使用 alglib::polynomialfit使用二阶多项式( m=3 在 alglib 中):

    alglib::barycentricinterpolant p;
    alglib::polynomialfitreport rep2;
    for (size_t i = 0; i < numPoints; i++)
    {
        // substitution: new t is equal to old t^2
        t_array(i) = std::pow(t_array(i), 2);
    }
    alglib::polynomialfit(t_array, y_array, 3, info, p, rep2);
    alglib::real_1d_array powerCoeffs;
    alglib::polynomialbar2pow(p, powerCoeffs);
    std::cout << "Success: " << info << std::endl;
    std::cout << "Results: " << powerCoeffs(0) << ", " << powerCoeffs(1) << ", " << powerCoeffs(2) << std::endl;
    std::cout << "RMS: " << rep2.rmserror << std::endl;

示例结果(在这种情况下它们是相同的，但通常该算法表现更好):

Success: 1
Results: -1000.04, -5.45533, 1.19111
RMS: 0.479453

alglib::polynomialfit很好，因为它使用重心插值，这在数值上更稳定。