Dopasowanie modelu wielomianowego do danych w R.

Przeczytałem odpowiedzi na to pytanie i są one bardzo pomocne, ale potrzebuję pomocy szczególnie w R.

Mam przykładowy zestaw danych w R w następujący sposób:

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)

Chcę dopasować model do tych danych, żeby y = f(x). Chcę, żeby był to model wielomianowy trzeciego rzędu.

Jak mogę to zrobić w R?

Ponadto czy R może mi pomóc znaleźć najlepiej dopasowany model?

Aby uzyskać wielomian trzeciego rzędu w x (x ^ 3), możesz to zrobić

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

lub

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

Możesz dopasować wielomian dziesiątego rzędu i uzyskać prawie idealne dopasowanie, ale czy powinieneś?

EDYCJA: poly (x, 3) jest prawdopodobnie lepszym wyborem (patrz @hadley poniżej).

To, który model jest „najlepiej dopasowanym modelem”, zależy od tego, co rozumiesz przez „najlepszy”. R ma narzędzia, które mogą Ci pomóc, ale aby dokonać wyboru, musisz podać definicję „najlepszego”. Rozważmy następujące przykładowe dane i kod:

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')

Który z tych modeli jest najlepszy? można by argumentować dla każdego z nich (ale ja nie chciałbym używać fioletowego do interpolacji).

Jeśli chodzi o pytanie `` czy R może mi pomóc znaleźć najlepiej dopasowany model '', prawdopodobnie istnieje funkcja, która to zrobi, zakładając, że możesz określić zestaw modeli do przetestowania, ale byłoby to dobre pierwsze podejście dla zbioru n-1 wielomiany stopnia:

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

Uwagi

• Słuszność tego podejścia zależy od celów, założenia optimize()i AIC()jeśli AIC jest kryterium, które chcesz użyć,

• polyfit()może nie mieć jednego minimum. sprawdź to za pomocą czegoś takiego:

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

• Użyłem tej as.integer()funkcji, ponieważ nie jest dla mnie jasne, jak zinterpretować wielomian niecałkowity.

• aby przetestować dowolny zestaw równań matematycznych, rozważ program „Eureqa” omówiony przez Andrew Gelmana tutaj

Aktualizacja

Zobacz także stepAICfunkcję (w pakiecie MASS) do automatyzacji wyboru modelu.

Najłatwiejszym sposobem znalezienia najlepszego dopasowania w R jest zakodowanie modelu jako:

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

Po zastosowaniu regresji obniżania AIC

lm.s <- step(lm.1)