在 R 中拟合“围绕”数据点的曲线

Question

我有一个由一组点组成的数据集。这些点在平面上的分布方式使得它们可以大致以抛物线为界。我正在尝试找到一种将抛物线拟合到点边界的方法。

这是我目前拥有的：

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）是通过反复试验确定的。

代码绘制点以及拟合抛物线。

现在这个系统可以工作了……但只是在某些时候。有时它会产生完全错误的拟合，我猜在这些情况下，逻辑项的参数是不合适的。运行代码几次看看我的意思。

我确信必须有一种更强大的方法来解决这个问题。想法和建议？

.

Answer 1

我无法提供完整的答案。我唯一的临时想法是为优化算法提供更好的起点 - 希望您更接近您尝试优化的函数的局部最小值。

估算一个粗略的第一个版本相当简单。如果你把你的抛物线写成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 个案例，一切都很顺利。您的结果是否经过审核，还是必须自动正确运行？

---

编辑 2

我对如何更新目标函数以使其更可靠有一些想法。现在我没有时间制定一个完整的解决方案，但也许这些想法可以帮助你：

我们的日期在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。但这是典型的：您需要在两个不同的目标（拟合数据、最大化抛物线）之间进行折衷。哪一项更重要的重量必须由您提交。

Answer 2

进一步感谢 thilo 非常有用的建议并纠正了我的幼稚想法。根据 thilo 的建议，使用抛物线下的面积和合适的惩罚函数，下面的解决方案似乎有效。我也改用 L-BFGS-B 优化，因为它在小 N 的情况下表现更好。

parabola.objective <- function(p) {
    d = p[2] * (data$x - p[1])^2 + p[3] - data$y
    #
    area <- function(x) {
        p[2] / 3 * (x - p[1])^3 + p[3] * x
    }
    #
    sum(- area(max(data$x)) + area(min(data$x)) + 100 * ifelse(d > 0, d^2, 0))
}

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

# the key to getting this working with a small number of points is the
# optimisation method: BFGS works well with around 300 points or more
# but L-BFGS-B seems to perform better down to around 100 points.
#
O = optim(c(A, B, C), parabola.objective, method="L-BFGS-B")

par = O$par

A = par[1]
B = par[2]
C = par[3]

curve(parabola, add = TRUE, lty = "dashed")