R 的poly() 函数为数据拟合生成正交多项式。但是，我想使用 R 之外的回归结果（比如在 C++ 中），并且似乎没有办法获得每个正交多项式的系数。
注1：我不是指回归系数，而是正交多项式本身的系数），例如
中的p_i(x)
y = a0 + a1*p_1(x) + a2*p_2(x) + ...
注意 2：我知道 poly(x, n, raw=T) 强制 poly 返回非正交多项式，但我想回归正交多项式，所以这就是我要寻找的。p>
【问题讨论】：
标签： r
多项式是使用您创建的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++ 中的递归公式计算多项式。
【讨论】：
注1：我不是指回归系数，而是正交多项式本身的系数），例如
p_i(x)在y = a0 + a1*p_1(x) + a2*p_2(x) + ...
作为josliber mentions，函数是递归构造的，表示它们的最简洁方式是使用 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);
}
如果需要，我们可以轻松摆脱对犰狳的依赖。我在下面验证我们得到与 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
【讨论】：