R：使用 qr 分解将具有正交多项式的模型转换为函数

Posted 2023-03-04

【中文标题】R：使用 qr 分解将具有正交多项式的模型转换为函数

【英文标题】：R: Translate a model having orthogonal polynomials to a function using qr decomposition

【发布时间】：2015-10-06 01:51:21

【问题描述】：

我正在使用 R 创建具有正交多项式的线性回归模型。我的模型是：

fit=lm(log(UFB2_BITRATE_REF3) ~ poly(QPB2_REF3,2)  + B2DBSA_REF3,data=UFB) 

UFB2_FPS_REF1= 29.98 27.65 26.30 25.69 24.68 23.07 22.96 22.16 21.51 20.75 20.75 26.15 24.59 22.91 21.02 19.59 18.80 18.21 17.07 16.74 15.98 15.80
QPB2_REF1 = 36 34 32 30 28 26 24 22 20 18 16 36 34 32 30 28 26 24 22 20 18 16
B2DBSA_REF1 = DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DOFFSOFF DONSON   DONSON  DONSON   DONSON   DONSON   DONSON   DONSON   DONSON   DONSON   DONSON   DONSON
Levels: DOFFSOFF DONSON

对应的总结是：

Call:
lm(formula = log(UFB2_BITRATE_REF3) ~ poly(QPB2_REF3, 2) + B2DBSA_REF3, data = UFB)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0150795 -0.0058792  0.0006155  0.0049245  0.0120587 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   9.630e+00  3.302e-02  291.62  < 2e-16 ***
poly(QPB2_REF3, 2, raw = T)1 -4.385e-02  2.640e-03  -16.61 2.31e-12 ***
poly(QPB2_REF3, 2, raw = T)2 -1.827e-03  5.047e-05  -36.20  < 2e-16 ***
B2DBSA_REF3DONSON            -3.746e-02  3.566e-03  -10.51 4.16e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.008363 on 18 degrees of freedom
Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 
F-statistic: 8.134e+04 on 3 and 18 DF,  p-value: < 2.2e-16 

接下来，我想为这个模型创建一个函数 f(x)=a + bx + cx^2 + ....。我想在 R 中使用 Gram Schmidt 算法使用 qr 分解。

你有什么想法吗？提前谢谢！

【问题讨论】：

您的示例中没有正交多项式。另外，请提供一个可重现的示例。

@Roland 函数poly() 返回一个矩阵，其列是正交多项式的基。

如果您指定 raw = TRUE 作为输出（但不是代码）显示，则不会。

@Roland 我不想指定raw = TRUE。

有什么不想要的理由吗？

【参考方案1】：

我忽略了“我想在 R 中使用 Gram Schmidt 算法使用 qr 分解”，但要注意 poly() 使用 qr() 来计算其正交多项式。

我读到这个问题是想采用正交多项式poly(QPB2_REF3, 2, raw = FALSE) 的系数模型，并用QPB2_REF3 的幂代数表达它。这意味着将正交多项式poly(QPB2_REF3, 2, raw = FALSE)1、poly(QPB2_REF3, 2, raw = FALSE)2 通常表示为QPB2_REF3 的幂系数，而不是attr(, "coefs") 对象的attr(, "coefs") 中的“中心化和归一化常数”。

多年来，在各种 R 论坛中，其他人也提出了类似的要求，被告知可以： (a) 使用 poly.predict() 计算多项式，因此不需要传统的形式系数； (b) 参见代码和/或 Kennedy & Gentle (1980, pp. 343-4) 中的算法。

(a) 没有满足我的教学需求。 在 (b) 上，我可以看到如何计算给定 x 的多项式值，但我只是迷失在试图推断传统形式系数的代数中：-

Kennedy & Gentle 提到了“在 p(x) 中求解 x”，在我看来，这是 lm 的简单建议，并导致了下面 orth2raw() 中实施的真正可怕的方法。我完全接受必须有一种更好、更直接的方法来从居中和归一化常数中推导出常规形式系数，但我无法解决，而且这种方法似乎有效。

orth2raw <- function(x)
# x <- poly(.., raw=FALSE) has a "coefs" attribute "which contains 
# the centering and normalization constants used in constructing 
# the orthogonal polynomials". orth2raw returns the coefficents of
# those polynomials in the conventional form
#    b0.x^0 + b1.x^1 + b2.x^2 + ...
# It handles the coefs list returned by my modifications of 
# poly and polym to handle multivariate predictions  
    o2r <- function(coefs)
       Xmean <- coefs$alpha[1]
       Xsd <- sqrt(coefs$norm2[3]/coefs$norm2[2])
       X <- seq(Xmean-3*Xsd, Xmean+3*Xsd, length.out=degree+1)
       Y <- poly(X, degree = degree, coefs=coefs)
       Rcoefs <- matrix(0,degree, degree+1)
       for (i in 1:degree) Rcoefs[i,1:(i+1)] <- coef(lm(Y[,i] ~ poly(X, i, raw=TRUE) ))
       dimnames(Rcoefs) <- list(paste0("poly(x)", 1:degree), paste0("x^",0:degree))
       Rcoefs
      
   degree <- max(attr(x, "degree"))
   coefs <- attr(x, "coefs")
   if(is.list(coefs[[1]])) lapply(coefs, o2r) else o2r(coefs)