r - 将多项式模型拟合到 R 中的数据

Question

我已经阅读了这个问题的答案，它们很有帮助，但我需要帮助。

我在 R 中有一个示例数据集，如下所示：

x <- c(32,64,96,118,126,144,152.5,158)  
y <- c(99.5,104.8,108.5,100,86,64,35.3,15)

我想为这些数据拟合一个模型，以便y = f(x). 我希望它是一个三阶多项式模型。

我怎么能在 R 中做到这一点？

此外，R 可以帮助我找到最合适的模型吗？

Answer 1

要获得 x (x^3) 中的三阶多项式，您可以执行

lm(y ~ x + I(x^2) + I(x^3))

或者

lm(y ~ poly(x, 3, raw=TRUE))

您可以拟合 10 阶多项式并获得近乎完美的拟合，但您应该这样做吗？

编辑： poly(x, 3) 可能是一个更好的选择（见下面的@hadley）。

Answer 2

哪个模型是“最佳拟合模型”取决于您所说的“最佳”。R 有工具可以提供帮助，但您需要提供“最佳”的定义以在它们之间进行选择。考虑以下示例数据和代码：

x <- 1:10
y <- x + c(-0.5,0.5)

plot(x,y, xlim=c(0,11), ylim=c(-1,12))

fit1 <- lm( y~offset(x) -1 )
fit2 <- lm( y~x )
fit3 <- lm( y~poly(x,3) )
fit4 <- lm( y~poly(x,9) )
library(splines)
fit5 <- lm( y~ns(x, 3) )
fit6 <- lm( y~ns(x, 9) )

fit7 <- lm( y ~ x + cos(x*pi) )

xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1, data.frame(x=xx)), col='blue')
lines(xx, predict(fit2, data.frame(x=xx)), col='green')
lines(xx, predict(fit3, data.frame(x=xx)), col='red')
lines(xx, predict(fit4, data.frame(x=xx)), col='purple')
lines(xx, predict(fit5, data.frame(x=xx)), col='orange')
lines(xx, predict(fit6, data.frame(x=xx)), col='grey')
lines(xx, predict(fit7, data.frame(x=xx)), col='black')

这些模型中哪一个是最好的？可以为它们中的任何一个进行参数（但我不想使用紫色的插值）。

Answer 3

关于“R 能帮我找到最佳拟合模型”的问题，可能有一个函数可以做到这一点，假设您可以说明要测试的模型集，但这对于 n-1 集来说是一个很好的第一种方法多项式：

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)

笔记

• 这种方法的有效性将取决于您的目标、假设optimize()以及AIC()AIC 是否是您想要使用的标准，

• polyfit()可能没有一个最小值。用类似的东西检查这个：

  for (i in 2:length(x)-1) print(polyfit(i))
  

• 我使用该as.integer()函数是因为我不清楚如何解释非整数多项式。

• 为了测试任意一组数学方程，请考虑Andrew Gelman在这里审查的“Eureqa”程序

更新

另请参阅stepAIC自动选择模型的功能（在 MASS 包中）。

Answer 4

在 R 中找到最佳拟合的最简单方法是将模型编码为：

lm.1 <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + ...)

使用降压 AIC 回归后

lm.s <- step(lm.1)

Answer 5

例如，如果我们想拟合一个 2 次多项式，我们可以直接通过以下方式求解线性方程组来实现：

以下示例显示了如何使用上述方程拟合抛物线 y = ax^2 + bx + c 并将其与lm()多项式回归解进行比较。希望这将有助于某人的理解，

x <- c(32,64,96,118,126,144,152.5,158)  
y <- c(99.5,104.8,108.5,100,86,64,35.3,15)

x4 <- sum(x^4)
x3 <- sum(x^3)
x2 <- sum(x^2)
x1 <- sum(x)
yx1 <- sum(y*x)
yx2 <- sum(y*x^2)
y1 <- sum(y)

A <- matrix(c(x4, x3, x2,
              x3, x2, x1,
              x2, x1, length(x)), nrow=3, byrow=TRUE)
B <- c(yx2,
       yx1,
       y1)

coef <- solve(A, B)   # solve the linear system of equations, assuming A is not singular
coef1 <- lm(y ~ x + I(x^2))$coef  # solution with lm

coef
# [1] -0.01345808  2.01570523 42.51491582
rev(coef1)
#     I(x^2)       x         (Intercept) 
#     -0.01345808  2.01570523 42.51491582 

plot(x, y, xlim=c(min(x), max(x)), ylim=c(min(y), max(y)+10), pch=19)
xx <- seq(min(x), max(x), 0.01)
lines(xx, coef[1]*xx^2+coef[2]*xx+coef[3], col='red', lwd=3, lty=5)
lines(xx,  coef1[3]*xx^2+ coef1[2]*xx+ coef1[1], col='blue')
legend('topright', legend=c("solve", "lm"),
       col=c("red", "blue"), lty=c(5,1), lwd=c(3,1), cex=0.8,
       title="quadratic fit", text.font=4)