r создать функцию из коэффициентов lm - PullRequest
1 голос
/ 10 ноября 2019

Код ниже отображает проблему. Я выполнил полиномиальную регрессию степени четыре lm на фрейме данных. df, чтобы получить model4. Затем я создаю регрессионную функцию fhat4. Это работает как задумано.

Я хочу обобщить это до любой степени полинома. Итак, я использую poly для создания modeln. Это соответствует model4. Но я не могу создать соответствующую функцию fhatn. Возможно, это как-то связано с циклом for?

df <- structure(list(x = c(0.3543937637005, 0.674911001464352, 0.21966037643142, 
0.14723521983251, 0.36166316177696, 0.975983075099066, 0.539355604210868, 
0.294046462047845, 0.853777077747509, 0.634912414476275), y = c(0.0120776002295315, 
0.655085238162428, 0.310665819328278, 0.525274415733293, 0.938241509487852, 
0.520828885724768, 0.241615766659379, 0.724816955626011, 0.808277940144762, 
0.358921303786337)), .Names = c("x", "y"), row.names = c(NA, 
-10L), class = "data.frame")

############################################# 
model4 <- lm(y~x+I(x^2)+I(x^3)+I(x^4), data=df)

fhat4 <- function (x) {
  model4$coefficients[1]+
  model4$coefficients[2]*x+
  model4$coefficients[3]*x^2+
  model4$coefficients[4]*x^3+
  model4$coefficients[5]*x^4
  }

fhat4(2)

############################################# 
modeln <- lm(y~poly(x,4,raw=TRUE), data=df)

fhatn <- function (x) {
  fn <- 0
  for (i in 0:5){
    fn <- fn + modeln$coefficients[i+1]*x^i
  }
}

fhatn(4)

1 Ответ

2 голосов
/ 10 ноября 2019

ваш for цикл должен идти только от 0 до 4, а не до 5. Также ваша функция ничего не возвращает, так что вы можете добавить return(fn) в конце.

В любом случае, вы можете реализоватьта же функция без каких-либо циклов:

modeln <- lm(y ~ poly(x, 4, raw = TRUE), data = df)

fhatn <- function (x) {
  sum(x^(seq_along(coef(modeln)) - 1) * coef(modeln))
}

fhatn(2)
[1] -150.6643

Обратите внимание, что coef(modeln) является альтернативой modeln$coefficients.

Или, как сказал Винсент в комментариях, вы можете использовать функцию предсказания:

predict(modeln, newdata = data.frame(x = 2))
-150.6643 
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...