如何使用`poly()`生成正交多项式？如何理解返回的“coefs”？

Question

如何使用`poly()`生成正交多项式？如何理解返回的“coefs”？

rmatrixregressionlinear-regressionlm

13

我对正交多项式的理解是，它们的形式为：

y(x) = a1 + a2(x - c1) + a3(x - c2)(x - c3) + a4(x - c4)(x - c5)(x - c6)... 直到所需项数为止

其中a1、a2等为每个正交项的系数（在不同的拟合中会有所变化），而c1、c2等则是正交项内的系数，通过确定这些系数使得各项保持正交性（在使用相同的x值进行拟合时保持一致）

我了解到poly()函数可用于拟合正交多项式。以下是一个例子：

x = c(1.160, 1.143, 1.126, 1.109, 1.079, 1.053, 1.040, 1.027, 1.015, 1.004, 0.994, 0.985, 0.977) # abscissae not equally spaced

y = c(1.217395, 1.604360, 2.834947, 4.585687, 8.770932, 9.996260, 9.264800, 9.155079, 7.949278, 7.317690, 6.377519, 6.409620, 6.643426)

# construct the orthogonal polynomial
orth_poly <- poly(x, degree = 5)

# fit y to orthogonal polynomial
model <- lm(y ~ orth_poly)

我想提取系数a1、a2等，以及正交系数c1、c2等。我不确定如何做到这一点。我的猜测是

model$coefficients

返回第一组系数，但我在如何提取其他系数方面遇到了困难。也许可以在内部进行操作。

attributes(orth_poly)$coefs

非常感谢。

- pyg

1个回答

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

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

- Zheyuan Li · Accepted Answer

我刚刚意识到2年前有一个非常相关的问题：从R的poly()函数中提取正交多项式系数？。那里的答案仅仅解释了predict.poly的作用，但是我的答案给出了完整的图片。

第一部分： poly 如何表示正交多项式

我理解的正交多项式具有如下形式：

y(x) = a1 + a2(x - c1) + a3(x - c2)(x - c3) + a4(x - c4)(x - c5)(x - c6)... 直到所需项数为止

不，没有这样干净的形式。poly()生成的是首项系数为一的正交多项式，可以通过以下递归算法表示。这就是predict.poly生成线性预测矩阵的方式。令人惊讶的是，poly本身并不使用这样的递归，而是使用暴力方法：对普通多项式的模型矩阵进行QR分解以获得正交空间。然而，这等价于递归。

第二部分：解释 poly() 的输出

让我们来看一个例子。以您帖子中的x为例：

X <- poly(x, degree = 5)

#                 1           2           3            4           5
# [1,]  0.484259711  0.48436462  0.48074040  0.351250507  0.25411350
# [2,]  0.406027697  0.20038942 -0.06236564 -0.303377083 -0.46801416
# [3,]  0.327795682 -0.02660187 -0.34049024 -0.338222850 -0.11788140
# ...           ...          ...        ...          ...         ...
#[12,] -0.321069852  0.28705108 -0.15397819 -0.006975615  0.16978124
#[13,] -0.357884918  0.42236400 -0.40180712  0.398738364 -0.34115435
#attr(,"coefs")
#attr(,"coefs")$alpha
#[1] 1.054769 1.078794 1.063917 1.075700 1.063079
# 
#attr(,"coefs")$norm2
#[1] 1.000000e+00 1.300000e+01 4.722031e-02 1.028848e-04 2.550358e-07
#[6] 5.567156e-10 1.156628e-12

以下是各个属性的含义：

alpha[1] 表示 x 的均值 x_bar，即中心位置；
alpha - alpha[1] 得到 alpha0, alpha1, ..., alpha4（alpha5 在 poly 返回矩阵 X 之前已经被计算出来，但由于不会在 predict.poly 中使用而被舍弃）；
norm2 的第一个值始终为1，其余部分为 l0, l1, ..., l5，它们表示矩阵 X 的列平方和；其中的 l0 是被舍弃的项 P0(x - x_bar) 的列平方和，总是等于 n（即 x 的长度），而第一个 1 只是为了使递归在 predict.poly 中继续进行；
beta0, beta1, beta2, ..., beta_5 没有被返回，但可以通过 norm2[-1] / norm2[-length(norm2)] 来计算。

第三部分：使用 QR 分解和递归算法实现 poly

正如之前提到的，poly 不使用递归，而 predict.poly 则使用了递归。个人认为这种设计不太一致，不太容易理解其逻辑和原因。在此我提供一个名为 my_poly 的函数，如果 QR = FALSE，则使用递归生成矩阵；当 QR = TRUE 时，则是类似但不完全相同的实现方式。代码有详细的注释，有助于您理解两种方法。

## return a model matrix for data `x`
my_poly <- function (x, degree = 1, QR = TRUE) {
  ## check feasibility
  if (length(unique(x)) < degree)
    stop("insufficient unique data points for specified degree!")
  ## centring covariates (so that `x` is orthogonal to intercept)
  centre <- mean(x)
  x <- x - centre
  if (QR) {
    ## QR factorization of design matrix of ordinary polynomial
    QR <- qr(outer(x, 0:degree, "^"))
    ## X <- qr.Q(QR) * rep(diag(QR$qr), each = length(x))
    ## i.e., column rescaling of Q factor by `diag(R)`
    ## also drop the intercept
    X <- qr.qy(QR, diag(diag(QR$qr), length(x), degree + 1))[, -1, drop = FALSE]
    ## now columns of `X` are orthorgonal to each other
    ## i.e., `crossprod(X)` is diagonal
    X2 <- X * X
    norm2 <- colSums(X * X)    ## squared L2 norm
    alpha <- drop(crossprod(X2, x)) / norm2
    beta <- norm2 / (c(length(x), norm2[-degree]))
    colnames(X) <- 1:degree
    } 
  else {
    beta <- alpha <- norm2 <- numeric(degree)
    ## repeat first polynomial `x` on all columns to initialize design matrix X
    X <- matrix(x, nrow = length(x), ncol = degree, dimnames = list(NULL, 1:degree))
    ## compute alpha[1] and beta[1]
    norm2[1] <- new_norm <- drop(crossprod(x))
    alpha[1] <- sum(x ^ 3) / new_norm
    beta[1] <- new_norm / length(x)
    if (degree > 1L) {
      old_norm <- new_norm
      ## second polynomial
      X[, 2] <- Xi <- (x - alpha[1]) * X[, 1] - beta[1]
      norm2[2] <- new_norm <- drop(crossprod(Xi))
      alpha[2] <- drop(crossprod(Xi * Xi, x)) / new_norm
      beta[2] <- new_norm / old_norm
      old_norm <- new_norm
      ## further polynomials obtained from recursion
      i <- 3
      while (i <= degree) {
        X[, i] <- Xi <- (x - alpha[i - 1]) * X[, i - 1] - beta[i - 1] * X[, i - 2]
        norm2[i] <- new_norm <- drop(crossprod(Xi))
        alpha[i] <- drop(crossprod(Xi * Xi, x)) / new_norm
        beta[i] <- new_norm / old_norm
        old_norm <- new_norm
        i <- i + 1
        }
      }
    }
  ## column rescaling so that `crossprod(X)` is an identity matrix
  scale <- sqrt(norm2)
  X <- X * rep(1 / scale, each = length(x))
  ## add attributes and return
  attr(X, "coefs") <- list(centre = centre, scale = scale, alpha = alpha[-degree], beta = beta[-degree])
  X
  }

第四部分：解释my_poly的输出

X <- my_poly(x, 5, FALSE)

因此，生成的矩阵与poly生成的矩阵相同，因此不予考虑。但是属性不同。

#attr(,"coefs")
#attr(,"coefs")$centre
#[1] 1.054769

#attr(,"coefs")$scale
#[1] 2.173023e-01 1.014321e-02 5.050106e-04 2.359482e-05 1.075466e-06

#attr(,"coefs")$alpha
#[1] 0.024025005 0.009147498 0.020930616 0.008309835

#attr(,"coefs")$beta
#[1] 0.003632331 0.002178825 0.002478848 0.002182892

my_poly返回更明显的构造信息：

centre给出 x_bar = mean(x)；
scale给出列范数（通过poly返回的norm2的平方根）；
alpha给出alpha1，alpha2，alpha3，alpha4；
beta给出beta1，beta2，beta3，beta4。

第五节：用于my_poly的预测例程

由于 my_poly 返回不同的属性，因此 stats:::predict.poly 与 my_poly 不兼容。这是适当的例程 my_predict_poly：

## return a linear predictor matrix, given a model matrix `X` and new data `x`
my_predict_poly <- function (X, x) {
  ## extract construction info
  coefs <- attr(X, "coefs")
  centre <- coefs$centre
  alpha <- coefs$alpha
  beta <- coefs$beta
  degree <- ncol(X)
  ## centring `x`
  x <- x - coefs$centre
  ## repeat first polynomial `x` on all columns to initialize design matrix X
  X <- matrix(x, length(x), degree, dimnames = list(NULL, 1:degree))
  if (degree > 1L) {
    ## second polynomial
    X[, 2] <- (x - alpha[1]) * X[, 1] - beta[1]
    ## further polynomials obtained from recursion
    i <- 3
    while (i <= degree) {
      X[, i] <- (x - alpha[i - 1]) * X[, i - 1] - beta[i - 1] * X[, i - 2]
      i <- i + 1
      }
    }
  ## column rescaling so that `crossprod(X)` is an identity matrix
  X * rep(1 / coefs$scale, each = length(x))
  }

考虑一个例子：

set.seed(0); x1 <- runif(5, min(x), max(x))

和

stats:::predict.poly(poly(x, 5), x1)
my_predict_poly(my_poly(x, 5, FALSE), x1)

给出完全相同的结果预测矩阵：

#               1          2           3          4          5
#[1,]  0.39726381  0.1721267 -0.10562568 -0.3312680 -0.4587345
#[2,] -0.13428822 -0.2050351  0.28374304 -0.0858400 -0.2202396
#[3,] -0.04450277 -0.3259792  0.16493099  0.2393501 -0.2634766
#[4,]  0.12454047 -0.3499992 -0.24270235  0.3411163  0.3891214
#[5,]  0.40695739  0.2034296 -0.05758283 -0.2999763 -0.4682834

请注意，预测程序只采用现有的构造信息，而不重构多项式。

第6节：将poly和predict.poly视为黑盒子即可

通常情况下，不需要理解其中的所有细节。对于统计建模来说，了解poly构建模型拟合的多项式基础并在lmObject$coefficients中找到其系数即可。进行预测时，用户无需调用predict.poly，因为predict.lm会自动调用它。因此，可以将poly和predict.poly视为黑盒子，这完全没有问题。