Regresja wielomianowa w R - z dodatkowymi ograniczeniami na krzywej

Question

Wiem, jak wykonać podstawową regresję wielomianową w R. Jednakże, mogę użyć tylko nls lub , aby dopasować linię, która minimalizuje błąd z punktami.

Działa to przez większość czasu, ale czasami, gdy w danych występują luki pomiarowe, model staje się bardzo sprzeczny z intuicją. Czy istnieje sposób na dodanie dodatkowych ograniczeń?

Powtarzalne Przykład:

Chcę dopasować model do następujących złożonych danych (podobnie do moich danych rzeczywistych):

x <- c(0, 6, 21, 41, 49, 63, 166) 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 

Najpierw wykreślić ją.

library(ggplot2) 
points <- ggplot(df, aes(x,y)) + geom_point(size=4, col='red') 
points 

Wygląda na to, czy mamy podłączone te punkty z linii, to zmiana kierunku 3 razy, więc spróbujmy zamontowanie Quartic do niego.

lm <- lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4)) 
quartic <- function(x) lm$coefficients[5]*x^4 + lm$coefficients[4]*x^3 + lm$coefficients[3]*x^2 + lm$coefficients[2]*x + lm$coefficients[1] 

points + stat_function(fun=quartic) 

Wygląda jak model pasuje punkty całkiem dobrze ... z wyjątkiem, ponieważ nasze dane miał dużą lukę pomiędzy 63 i 166, istnieje ogromny skok tam, który ma powody do w modelu. (Dla mojego danych rzeczywistych wiem, że nie ma tam ogromny szczyt)

Więc pytanie jest w tym przypadku:

• Jak mogę ustawić, że lokalne maksimum być na (166, 9,8)?

Jeśli nie jest to możliwe, to kolejny sposób, aby zrobić to byłoby:

• Jak mogę ograniczyć Y wartości przewidywanych przez linię z coraz większym niż y = 9,8.

A może jest lepszy model do użycia? (Poza tym, że robi to kawałek po kawałku). Moim celem jest porównanie cech modeli między wykresami.

Answer 1

Typ funkcji doskonale pasuje do twoich danych (ale nie dla celów predykcji). Krzywe splajnu są szeroko stosowane w obszarach CAD i czasami po prostu pasują do danych w matematyce i mogą oznaczać brak znaczenia fizyki w porównaniu z regresją. Więcej informacji pod numerem here i wprowadzenie tła w here.

example(spline) pokaże Ci mnóstwo fantazyjnych przykładów, a ja faktycznie używam jednego z nich.

Ponadto, będzie bardziej uzasadnione próbki więcej punktów danych, a następnie dopasowanie przez lm lub nls regresji dla predykcji.

Przykładowy kod:

library(splines) 

x <- c(0, 6, 21, 41, 49, 63, 166) 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 

s1 <- splinefun(x, y, method = "monoH.FC") 

plot(x, y) 
curve(s1(x), add = TRUE, col = "red", n = 1001) 

Innym podejściem mogę myśl jest przymusu zakres parametrów w regresji, dzięki czemu można uzyskać przewidywane dane w swoim oczekiwanym zakresie.

Bardzo prosty kod z optim w poniżej, ale tylko wybór.

dat <- as.data.frame(cbind(x,y)) 
names(dat) <- c("x", "y") 

# your lm 
# lm<-lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4)) 

# define loss function, you can change to others 
min.OLS <- function(data, par) { 
     with(data, sum(( par[1]  + 
         par[2] * x + 
         par[3] * (x^2) + 
         par[4] * (x^3) + 
         par[5] * (x^4) + 
         - y)^2) 
      ) 
} 

# set upper & lower bound for your regression 
result.opt <- optim(par = c(0,0,0,0,0), 
       min.OLS, 
       data = dat, 
       lower=c(3.6,-2,-2,-2,-2), 
       upper=c(6,1,1,1,1), 
       method="L-BFGS-B" 
) 

predict.yy <- function(data, par) { 
       print(with(data, ((
        par[1]  + 
        par[2] * x + 
        par[3] * (x^2) + 
        par[4] * (x^3) + 
        par[5] * (x^4)))) 
       ) 
    } 

    plot(x, y, main="LM with constrains") 
    lines(x, predict.yy(dat, result.opt$par), col="red") 

Answer 2

pójdę do regresji lokalnym jak eipi10 sugerowane. Jednakże, jeśli chcesz uzyskać regresję wielomianu, możesz spróbować zminimalizować karę sumy kwadratów.

Oto przykład, w którym funkcja popełnia odbiegające „zbyt dużo” od linii prostej:

library(ggplot2) 
library(maxLik) 
x <- c(0, 6, 21, 41, 49, 63, 166)/100 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 
points <- ggplot(df, aes(x,y)) + geom_point(size=4, col='red') 

polyf <- function(par, x=df$x) { 
    ## the polynomial function 
    par[1]*x + par[2]*x^2 + par[3]*x^3 + par[4]*x^4 + par[5] 
} 
quarticP <- function(x) { 
    polyf(par, x) 
} 
## a evenly distributed set of points, penalize deviations on these 
grid <- seq(range(df$x)[1], range(df$x)[2], length=10) 

objectiveF <- function(par, kappa=0) { 
    ## Calculate penalized sum of squares: penalty for deviating from linear 
    ## prediction 
    PSS <- sum((df$y - polyf(par))^2) + kappa*(pred1 - polyf(par))^2 
    -PSS 
} 

## first compute linear model prediction 
res1 <- lm(y~x, data=df) 
pred1 <- predict(res1, newdata=data.frame(x=grid)) 
points <- points + geom_smooth(method='lm',formula=y~x) 
print(points) 

## non-penalized function 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0)) 
par <- coef(res) 
points <- points + stat_function(fun=quarticP, col="green") 
print(points) 

## penalty 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0), kappa=0.5) 
par <- coef(res) 
points <- points + stat_function(fun=quarticP, col="yellow") 
print(points) 

Wynik z karą 0,5 wyglądem wygląda następująco: Można dostosować karę, a grid, miejsca, w których odstępstwa są karane.

Answer 3

Źródło Ott Toomets nie działa dla mnie, były pewne błędy. Oto wersja poprawiona (bez użycia ggplot2):

library(maxLik) 
x <- c(0, 6, 21, 41, 49, 63, 166)/100 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 

polyf <- function(par, x=df$x) { 
    ## the polynomial function 
    par[1]*x + par[2]*x^2 + par[3]*x^3 + par[4]*x^4 + par[5] 
} 
quarticP <- function(x) { 
    polyf(par, x) 
} 
## a evenly distributed set of points, penalize deviations on these 
grid <- seq(range(df$x)[1], range(df$x)[2], length=10) 

objectiveF <- function(par, kappa=0) { 
    ## Calculate penalized sum of squares: penalty for deviating from linear 
    ## prediction 
    PSS <- sum((df$y - polyf(par))^2) + kappa*(pred1 - polyf(par, x=grid))^2 
    -PSS 
} 

plot(x,y, ylim=c(0,10)) 

## first compute linear model prediction 
res1 <- lm(y~x, data=df) 
pred1 <- predict(res1, newdata=data.frame(x=grid)) 
coefs = coef(res1) 
names(coefs) = NULL 
constant = coefs[1] 
xCoefficient = coefs[2] 
par = c(xCoefficient,0,0,0,constant) 

curve(quarticP, from=0, to=2, col="black", add=T) 


## non-penalized function 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0)) 
par <- coef(res) 
curve(quarticP, from=0, to=2, col="red", add=T) 

## penalty 
res2 <- maxBFGS(objectiveF, start=c(0,0,0,0,0), kappa=0.5) 
par <- coef(res2) 
curve(quarticP, from=0, to=2, col="green", add=T)