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

Question

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返回非正交多项式，但是我想对正交多项式进行回归，所以这就是我要寻找的。

Answer 1

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

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

为了表示，让我们将a_d的索引d中的元素称为alpha，并将我们将n_d的索引d中的元素称为norm2。 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 ++中使用递归公式来评估多项式。