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

Question

R poly()函数产生用于数据拟合的正交多项式.但是,我想使用R之外的回归结果(比如在C++中),并且似乎没有办法获得每个正交多项式的系数.

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

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

注2:我知道poly(x, n, raw=T)强制多边形来返回非正交多项式,但我想对正交多项式进行回归,这就是我正在寻找的东西.

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++来评估你的多项式.