Полиномиальная регрессия в R - с дополнительными ограничениями на кривую

Я знаю, как сделать базовую полиномиальную регрессию в R. Однако я могу использовать nls или lm для соответствия строке, которая минимизирует ошибку с точками.

Это работает большую часть времени, но иногда, когда в данных имеются пробелы измерения, модель становится очень противоречивой. Есть ли способ добавить дополнительные ограничения?

Воспроизводимый пример:

Я хочу подгонять модель к следующим составленным данным (подобно моим реальным данным):

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)

Сначала запишите его.

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

Похоже, что если бы мы связали эти точки с линией, это изменило бы направление 3 раза, поэтому попробуйте установить для него четверть.

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)

Похоже, что модель очень хорошо подходит к точкам... кроме того, потому что у наших данных был большой разрыв между 63 и 166, там есть огромный шип, который не имеет причины находиться в модели. (Для моих фактических данных я знаю, что там нет огромного пика)

Итак, вопрос в этом случае:

Как установить максимальный локальный максимум (166, 9.8)?

Если это невозможно, то другой способ сделать это:

Как ограничить значения y, предсказанные линией, больше, чем y = 9.8.

Или, может быть, есть лучшая модель для использования? (Кроме того, что делать это кусочно). Моя цель - сравнить особенности моделей между графиками.

Ответ 1

Функция spline будет идеально соответствовать вашим данным (но не предназначена для целей прогнозирования). Кривые сплайна широко используются в областях САПР, и когда-то это просто соответствует точке данных в математике и может быть недостатком физического смысла по сравнению с регрессией. Дополнительная информация в здесь и отличное введение в здесь.

example(spline) покажет вам множество причудливых примеров, и на самом деле я использую один из них.

Кроме того,, будет более разумно отбирать больше точек данных, а затем подстраиваться регрессией lm или nls для прогнозирования.

Пример кода:

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)

Другим подходом, который я могу представить, является ограничение диапазона параметров в регрессии, чтобы вы могли получить предсказанные данные в ожидаемом диапазоне.

Очень простой код с optim ниже, но только выбор.

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

Ответ 2

Я бы пошел на местную регрессию, как предположил eipi10. Однако, если вы хотите иметь полиномиальную регрессию, вы можете попытаться свести к минимуму штрафную сумму квадратов.

Вот пример, когда функция штрафуется за отклонение "слишком много" от прямой:

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)

Результат с штрафом 0.5 выглядит следующим образом: Вы можете отрегулировать штраф, а grid - места, где отклонения будут оштрафованы.

Ответ 3

Источник Ott Toomets не работал у меня, были некоторые ошибки. Вот исправленная версия (без использования 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)