R: Преобразуйте модель, имеющую ортогональные полиномы, в функцию, используя разложение qr.

Я использую 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 + .... для этой модели. Я хочу использовать разложение qr с использованием алгоритма Грама Шмидта в R.

У тебя есть что-нибудь на уме? Заранее спасибо!


r polynomials orthogonal
person zinon    schedule 16.07.2015    source источник
comment
В вашем примере нет ортогональных многочленов. Кроме того, предоставьте воспроизводимый пример.   -  person Roland    schedule 16.07.2015
comment
@Roland Функция poly() возвращает матрицу, столбцы которой являются основой ортогональных многочленов.   -  person zinon    schedule 16.07.2015
comment
Нет, если вы укажете raw = TRUE, как показывает ваш вывод (но не код).   -  person Roland    schedule 16.07.2015
comment
@ Роланд, я не хочу указывать raw = TRUE.   -  person zinon    schedule 16.07.2015
comment
Есть причина, по которой ты этого не хочешь?   -  person Roland    schedule 16.07.2015
comment
@Roland Я провожу исследование и, используя ортогональные полиномы, получаю более точные прогнозы.   -  person zinon    schedule 16.07.2015
comment
Что ж, если вам нужна помощь, исправьте свой пример, чтобы он показывал правильный вывод и (что более важно) мог быть воспроизведен, просто скопировав и запустив его в новом сеансе R.   -  person Roland    schedule 16.07.2015
comment
@Roland Что вы подразумеваете под правильным выводом?   -  person zinon    schedule 16.07.2015
comment
Правильный вывод означает вывод, который действительно показывает результат с ортогональными полиномами. В любом случае, вы не получите функцию, подобную f(x)=a + bx + cx^2 + ..... Вместо этого вы получите f(x) = g(x) %*% coef(fit), где g(x) создает матрицу дизайна.   -  person Roland    schedule 16.07.2015
comment
Я не знаю, как выглядит результат, так, как я хочу его создать.   -  person zinon    schedule 16.07.2015
comment
@Roland Я добавил свои данные к своему вопросу.   -  person zinon    schedule 16.07.2015
comment
Добавьте его более полезным способом. Прочитайте stackoverflow.com/a/5963610/1412059.   -  person Roland    schedule 16.07.2015

Ответы (1)


arrow_upward
2
arrow_downward

Я игнорирую «Я хочу использовать разложение qr с использованием алгоритма Грама Шмидта в R», за исключением того, что 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") объекта poly().

На протяжении многих лет на различных форумах R другие люди делали аналогичные запросы, чтобы им сказали, что можно: (а) вычислить полиномы с помощью poly.predict(), поэтому коэффициенты обычной формы не нужны; (b) см. алгоритм в коде и/или Kennedy & Gentle (1980, стр. 343–344).

а) не отвечал моим дидактическим потребностям. На (b) я мог видеть, как вычислить значения полинома для заданного x, но я просто заблудился в алгебре, пытаясь вывести коэффициенты обычной формы :-{

Кеннеди и Джентл ссылаются на «нахождение x в p(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)
   }
person user20637    schedule 17.07.2015