Вы должны минимизировать сумму квадратов с учетом ограничения; lm не допускает подобных ограничений, поэтому вы должны использовать общую функцию оптимизации, такую как optim. Вот один из способов сделать это.
Составьте некоторые данные. Здесь я скажу, что известный максимум 50.
set.seed(5)
d <- data.frame(x=seq(-5, 5, len=51))
d$y <- 50 - 0.3*d$x^2 + rnorm(nrow(d))
M <- 50
Создать функцию для получения квадратичной кривой для точек в точке x
с заданными квадратичными и линейными коэффициентами и заданным максимумом М.
Исчисление является простым; подробности см. в ответе Даффимо.
qM <- function(a, b, x, M) {
c <- M - (3*b^2)/(4*a)
a*x^2 + b*x + c
}
Создайте функцию, которая получит сумму квадратов между
квадратичная кривая с заданными квадратичными и линейными коэффициентами
и данные в д.
ff <- function(ab, d, M) {
p <- qM(ab[1], ab[2], d$x, M)
y <- d$y
sum((p-y)^2)
}
Получите обычный lm, пригодный для использования в качестве начальных значений.
m0 <- lm(y ~ I(x^2) + x, data=d)
start <- coef(m0)[2:3]
Оптимизация коэффициентов в функции ff.
o <- optim(start, ff, d=d, M=M)
o$par
Составьте график, показывающий, как подгонка имеет максимум при 50; оригинальная lm не подходит.
plot(d)
xs <- seq(-5, 5, len=101)
lines(xs, predict(m0, newdata=data.frame(x=xs)), col="gray")
lines(xs, qM(o$par[1], o$par[2], xs, M))
abline(h=50, lty=3)
