我有一个由一组点组成的数据集.这些点以这样一种方式分布在飞机上,即它们可以用抛物线大致界定.我试图找到一种方法将抛物线拟合到点的边界.
这就是我目前所拥有的:
a = 1
b = 2
c = 3
parabola <- function(x) {
a * x^2 + b * x + c
}
N = 10000
x <- runif(N, -4, 3)
y <- runif(N, 0, 10)
data <- data.frame(x, y)
data <- subset(data, y >= parabola(x))
plot(data, xlim = c(-5, 5), ylim = c(0, 10), col = "grey")
fr <- function(x) {
PAR = x[1] * data$x^2 + x[2] * data$x + x[3]
#
sum((PAR - data$y)^2 + 100 * plogis(PAR - data$y, scale = 0.00001))
}
par = optim(c(0, 0, 0), fr)$par
a = par[1]
b = par[2]
c = par[3]
curve(parabola, add = TRUE, lty = "dashed")
这将创建一个样本数据集,然后将曲线拟合到边界.目标函数包括一个"正常"平方误差项,它适用于抛物线数据,以及第二个逻辑项,它惩罚生活在抛物线下方的点.第二项的参数(100和0.00001)通过反复试验确定.
代码绘制点以及拟合的抛物线.
现在这个系统工作......但只有一些时间.有时它产生一个完全错误的拟合,我想在这些情况下,逻辑术语的参数是不合适的.运行代码几次,看看我的意思.
我相信必须有一种更强大的方法来解决这个问题.想法和建议？
.
我无法提供完整的答案。我唯一的临时想法是为优化算法提供更好的起点 - 希望您更接近您尝试优化的函数的局部最小值。
估计粗略的第一个版本相当简单。b*(x-a)^2+c
如果你写出你可以估计的抛物线
a <- data$x[which.min(data$y)]
c <- min(data$y)
b1 <- (data$y[which.min(data$x)] - c) / (min(data$x) - a)^2
b2 <- (data$y[which.max(data$x)] - c) / (max(data$x) - a)^2
b <- mean(c(b1, b2))
我根据我的建议和“BFGS”方法进行了另一次密集测试。我找不到采用以下方法的反例：
seed <- floor(runif(1,1,1000))
set.seed(seed)
a = 1
b = 2
c = 3
parabola <- function(x) {
b * (x-a)^2 + c
}
N = 10000
x <- runif(N, -4, 3)
y <- runif(N, 0, 10)
data <- data.frame(x, y)
data <- subset(data, y >= parabola(x))
plot(data, xlim = c(-5, 5), ylim = c(0, 10), col = "grey")
fr <- function(x) {
PAR = x[2] * (data$x - x[1])^2 + x[3]
#
sum((PAR - data$y)^2 + 100 * plogis(PAR - data$y, scale = 0.00001))
}
a <- data$x[which.min(data$y)]
c <- min(data$y)
b1 <- (data$y[which.min(data$x)] - c) / (min(data$x) - a)^2
b2 <- (data$y[which.max(data$x)] - c) / (max(data$x) - a)^2
b <- mean(c(b1, b2))
par = optim(c(a, b, c), fr, method="BFGS")$par
a = par[1]
b = par[2]
c = par[3]
curve(parabola, add = TRUE, lty = "dashed")
但是，不能保证正确收敛。我尝试了大约50个案例，一切都很顺利。您的结果是否经过审核或者是否必须自动正确运行？
我对如何更新目标函数以使其更加可靠有一些想法。现在我没有时间制定完整的解决方案，但也许这个想法可以帮助你：
我们的日期在 内range(data$x)。现在我们想要找到一条尽可能适合该数据下边界的抛物线，或者换句话说，找到最大化的值 a、b、c
\int_{\range(x)} ax^2 + bx+c dx
（请原谅笨拙的 LaTeX - 有时编写公式更好）。
现在，可以使用像这样的惩罚函数来惩罚抛物线以下的点
\lambda (ax_i^2+bx_i+c - y_i)^2 if below parabola, 0 otherwise
从区间中减去该函数应该会得到一个合适的、平滑的目标函数。尽可能简化函数似乎是比使用最小二乘法更好的模型，最小二乘法试图拟合一条穿过数据点中间的线。
不过，您仍然需要选择合适的 lambda。但这是典型的：您需要在两个不同目标之间进行折衷（拟合数据、最大化抛物线）。哪一个更重要必须由您提交。