在R中“围绕”数据点拟合曲线
我有一个由点集合组成的数据集。这些点以这样的方式分布在平面上，即它们可以粗略地以抛物线为边界。我正在试着找到一种方法来拟合抛物线到点的边界。
以下是我目前的资料：
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")这将创建一个样本数据集，然后将曲线拟合到边界。目标函数由一个“正常”的平方误差项组成，它拟合数据的抛物线，以及第二个logistic项，它惩罚生活在抛物线以下的点。该第二项的参数(100和0.00001)是通过试错法确定的。
代码将绘制点以及拟合的抛物线。
现在这个系统可以工作了..。但只有在某些情况下。有时它会产生完全错误的拟合，我猜在这些情况下，逻辑术语的参数是不合适的。运行代码几次，以了解我的意思。
我相信一定有更强大的方法来解决这个问题。想法和建议？
。
回答 2
Stack Overflow用户
发布于 2012-11-23 15:40:31
我不能提供一个完整的答案。我唯一特别的想法是为优化算法提供更好的起点-希望你更接近你试图优化的函数的局部最小值。
估计一个粗略的第一个版本相当简单。如果你把抛物线写成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。但这是典型的:你需要在两个不同的目标之间达成妥协(拟合数据，最大化抛物线)。哪一个更重要的权重必须由你提交。
Stack Overflow用户
发布于 2012-11-26 15:03:58
进一步感谢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")https://stackoverflow.com/questions/13524080
复制相似问题
腾讯云开发者
