在R中将曲线“套”在数据点周围

Question

在R中将曲线“套”在数据点周围

8

我有一个数据集，其中包含一组点。这些点在平面上分布，可以大致被一个抛物线界定。我正在寻找一种方法来拟合这些点的边界所形成的抛物线。

目前我拥有以下内容：

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

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

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

我相信一定有更稳健的方法来解决这个问题。有什么想法和建议吗？

.

- datawookie

optim中的默认算法不是很好。尝试指定method="BFGS"或method="L-BFGS-B". - Hong Ooi

将随机数种子设置为可重现的问题示例。如果我执行set.seed(999)并运行您的代码，那是否是一个错误拟合的例子？给我们一些可处理的东西！在这种情况下，method="BFGS"会产生更好的拟合效果... - Spacedman

是的，确切地说：使用种子值为999会失败，但使用种子值为1时可以正常工作。 - datawookie

2个回答

0

特别感谢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")

- datawookie

我必须承认我有点惊讶。使用您的优化方法，拟合比给定点更多的抛物线应该产生相同的目标函数值。我能想象的唯一原因是起始值通常适合一个“太小”的抛物线，而优化恰好命中最佳拟合函数附近的函数。我仍然建议添加一些罚项来限制过大的抛物线（即使只是 +p[3] + c*p[2]，这可能已经足够了）。 - Thilo

嗨Thilo，你说得很对。这确实适用于由我的示例代码生成的所有漂亮，整洁的“测试”案例。但是当我将其带到现实世界中时，它却出了大问题。我太天真了，以为这会起作用。如果初始抛物线在数据点之间开始，则此目标函数似乎会产生合理的结果。否则，它会彻底失败。所以你的惊讶完全是有道理的，而我则感到谦卑。 - datawookie

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Thilo · Accepted Answer

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

估计一个粗略的第一个版本相当简单。如果你把你的抛物线写成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

从该函数中减去区间应该会给你一个合适、平滑的目标函数。尽可能简化函数似乎比使用最小二乘法更好，最小二乘法试图通过数据点的“中间”拟合一条线。

不过，你仍然需要选择一个合适的λ。但这是典型的：你需要在两个不同的目标之间做出妥协（拟合数据，最大化抛物线）。哪个更重要的权重必须由你提交。