从R的poly()函数中提取正交多项式系数？

Question

从R的poly()函数中提取正交多项式系数？

19

R的poly()函数可以用于数据拟合产生正交多项式。但是，我想在R之外（比如在C++中）使用回归的结果，目前似乎无法获得每个正交多项式的系数。

注1：我指的不是回归系数，而是正交多项式本身的系数，例如 p_i(x).

y = a0 + a1*p_1(x) + a2*p_2(x) + ...

注意2: 我知道poly(x, n, raw=T)会强制多项式返回非正交多项式，但我想回归正交多项式，因此这就是我要寻找的。

- Gilead

2个回答

2

Note 1: I don't mean the regression coefficients, but the coefficients of the orthogonal polynomials themselves), e.g. the p_i(x) in
y = a0 + a1*p_1(x) + a2*p_2(x) + ...

正如josliber提到的，这些函数是递归构建的，最紧凑的表示方法是使用R使用的规范和“alpha”系数。使用（Rcpp）Armadillo的C ++版本可能类似于：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

namespace polycpp {
class poly_basis {
  void scale_basis(arma::mat &X) const {
    if(X.n_cols > 0)
      X.each_row() /= 
        arma::sqrt(norm2.subvec(1L, norm2.n_elem - 1L)).t();
  }

public:
  arma::vec const alpha, norm2;

  poly_basis(arma::vec const &alpha, arma::vec const &norm2):
    alpha(alpha), norm2(norm2) { 
    for(size_t i = 0; i < norm2.size(); ++i)
      if(norm2[i] <= 0.)
        throw std::invalid_argument("invalid norm2");
    if(alpha.n_elem + 2L != norm2.n_elem)
      throw std::invalid_argument("invalid alpha");
  }

  /**
   behaves like poly(x, degree). The orthogonal polynomial is returned by 
   reference.
   */
  static poly_basis get_poly_basis
    (arma::vec x, size_t const degree, arma::mat &X){
    size_t const n = x.n_elem, 
                nc = degree + 1L;
    double const x_bar = arma::mean(x); 
    x -= x_bar;
    arma::mat XX(n, nc);
    XX.col(0).ones();
    for(size_t d = 1L; d < nc; d++){
      double       * xx_new = XX.colptr(d);
      double const * xx_old = XX.colptr(d - 1);
      for(size_t i = 0; i < n; ++i, ++xx_new, ++xx_old)
        *xx_new = *xx_old * x[i];
    }

    arma::mat R;
    /* TODO: can be done smarter by calling LAPACK or LINPACK directly */
    if(!arma::qr_econ(X, R, XX))
      throw std::runtime_error("QR decomposition failed");

    for(size_t c = 0; c < nc; ++c)
      X.col(c) *= R.at(c, c);

    arma::vec norm2(nc + 1L), 
              alpha(nc - 1L);
    norm2[0] = 1.;
    for(size_t c = 0; c < nc; ++c){
      double z_sq(0), 
           x_z_sq(0);
      double const *X_i = X.colptr(c);
      for(size_t i = 0; i < n; ++i, ++X_i){
        double const z_sq_i = *X_i * *X_i;
        z_sq += z_sq_i;
        if(c < degree)
          x_z_sq += x[i] * z_sq_i;
      }
      norm2[c + 1] = z_sq;
      if(c < degree)
        alpha[c] = x_z_sq / z_sq + x_bar;
    }

    poly_basis out(alpha, norm2);
    out.scale_basis(X);
    return out;
  }

  /** behaves like predict(<poly>, newdata). */
  arma::mat operator()(arma::vec const &x) const {
    size_t const n = x.n_elem;
    arma::mat out(n, alpha.n_elem + 1L);
    out.col(0).ones();
    if(alpha.n_elem > 0L){
      out.col(1) = x;
      out.col(1) -= alpha[0];
      for(size_t c = 1; c < alpha.n_elem; c++){
        double       * x_new  = out.colptr(c + 1L);
        double const * x_prev = out.colptr(c), 
                     * x_old  = out.colptr(c - 1L), 
                     * x_i    = x.memptr();
        double const fac = norm2[c + 1L] / norm2[c];
        for(size_t i = 0; i < n; ++i, ++x_new, ++x_prev, ++x_old, ++x_i)
          *x_new = (*x_i - alpha[c]) * *x_prev - fac * *x_old;
      } 
    }

    scale_basis(out);
    return out;
  }
};
} // namespace polycpp

// export the functions to R to show that we get the same
using namespace polycpp;
using namespace Rcpp;

// [[Rcpp::export(rng = false)]]
List my_poly(arma::vec const &x, unsigned const degree){
  arma::mat out;
  auto basis = poly_basis::get_poly_basis(x, degree, out);
  return List::create(
    Named("X") = out, 
    Named("norm2") = basis.norm2, 
    Named("alpha") = basis.alpha);
}

// [[Rcpp::export(rng = false)]]
arma::mat my_poly_predict(arma::vec const &x, arma::vec const &alpha, 
                          arma::vec const &norm2){
  poly_basis basis(alpha, norm2);
  return basis(x);
}

如果需要，我们可以轻松摆脱Armadillo依赖。我在下面验证，我们得到与R函数相同的结果。

set.seed(1L)
x <- rnorm(100)
Rp <- poly(x, degree = 4L)
Cp <- my_poly(x, 4L)

all.equal(unclass(Rp), Cp$X[, -1L], check.attributes = FALSE)
#R> [1] TRUE
all.equal(attr(Rp, "coefs"), 
          lapply(Cp[c("alpha", "norm2")], drop))
#R> [1] TRUE

z <- rnorm(20)
Rpred <- predict(Rp, z)
Cpred <- my_poly_predict(z, Cp$alpha, Cp$norm2)
all.equal(Rpred, Cpred[, -1], check.attributes = FALSE)
#R> [1] TRUE
all.equal(Cp$X, my_poly_predict(x, Cp$alpha, Cp$norm2))
#R> [1] TRUE

一个不错的奖励，尽管在实践中可能并不重要，是新功能更快。

options(digits = 3)
microbenchmark::microbenchmark(
  R = poly(x, degree = 4L), cpp = my_poly(x, 4L))
#R> Unit: microseconds
#R> expr    min     lq  mean median    uq   max neval
#R>    R 118.93 123.63 135.4  126.1 129.0 469.1   100
#R>  cpp   7.22   7.97  11.3   10.9  11.7  89.4   100
microbenchmark::microbenchmark(
  R = predict(Rp, z), cpp = my_poly_predict(z, Cp$alpha, Cp$norm2))
#R> Unit: microseconds
#R> expr  min    lq  mean median    uq   max neval
#R>    R 18.6 19.20 20.50  19.43 19.89 92.86   100
#R>  cpp  1.2  1.39  1.92   1.98  2.23  8.85   100

- Benjamin Christoffersen

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- josliber · Accepted Answer

多项式的定义是使用您创建的poly对象的alpha和norm2系数进行递归定义的。让我们看一个例子：

z <- poly(1:10, 3)
attributes(z)$coefs
# $alpha
# [1] 5.5 5.5 5.5
# $norm2
# [1]    1.0   10.0   82.5  528.0 3088.8

为了表示方便，我们将alpha中索引为d的元素称为a_d，将norm2中索引为d的元素称为n_d。 F_d(x)是生成的度数为d的正交多项式。对于一些基本情况，我们有：

F_0(x) = 1 / sqrt(n_2)
F_1(x) = (x-a_1) / sqrt(n_3)

其余的多项式是通过递归定义的：

F_d(x) = [(x-a_d) * sqrt(n_{d+1}) * F_{d-1}(x) - n_{d+1} / sqrt(n_d) * F_{d-2}(x)] / sqrt(n_{d+2})

确认 x=2.1:

x <- 2.1
predict(z, newdata=x)
#               1         2         3
# [1,] -0.3743277 0.1440493 0.1890351
# ...

a <- attributes(z)$coefs$alpha
n <- attributes(z)$coefs$norm2
f0 <- 1 / sqrt(n[2])
(f1 <- (x-a[1]) / sqrt(n[3]))
# [1] -0.3743277
(f2 <- ((x-a[2]) * sqrt(n[3]) * f1 - n[3] / sqrt(n[2]) * f0) / sqrt(n[4]))
# [1] 0.1440493
(f3 <- ((x-a[3]) * sqrt(n[4]) * f2 - n[4] / sqrt(n[3]) * f1) / sqrt(n[5]))
# [1] 0.1890351

将多项式导出到您的C++代码中最紧凑的方法可能是导出attributes(z)$coefs$alpha和attributes(z)$coefs$norm2，然后在C++中使用递归公式来评估您的多项式。