在曲线上找到最佳权衡点

Posted 2023-03-10

【中文标题】在曲线上找到最佳权衡点

【英文标题】：Finding the best trade-off point on a curve

【发布时间】：2011-01-02 08:31:56

【问题描述】：

假设我有一些数据，我想为其拟合一个参数化模型。我的目标是找到这个模型参数的最佳值。

我正在使用AIC/BIC/MDL 类型的标准进行模型选择，该标准奖励具有低错误的模型但也惩罚具有高复杂性的模型（我们正在寻找最简单但最有说服力的解释这个数据可以这么说，Occam's razor)。

根据上述情况，这是我根据三个不同标准（两个要最小化，一个要最大化）得到的东西的一个例子：

您可以很容易地看到肘部形状，并且可以在该区域的某处为参数选择一个值。 问题是我这样做是为了进行大量实验，我需要一种无需干预即可找到此值的方法。

我的第一个直觉是尝试从拐角处以 45 度角画一条线并继续移动它直到它与曲线相交，但这说起来容易做起来难 :) 如果曲线是有点歪。

关于如何实现这一点或更好的想法有什么想法吗？

这里是重现上述图之一所需的样本：

curve = [8.4663 8.3457 5.4507 5.3275 4.8305 4.7895 4.6889 4.6833 4.6819 4.6542 4.6501 4.6287 4.6162 4.585 4.5535 4.5134 4.474 4.4089 4.3797 4.3494 4.3268 4.3218 4.3206 4.3206 4.3203 4.2975 4.2864 4.2821 4.2544 4.2288 4.2281 4.2265 4.2226 4.2206 4.2146 4.2144 4.2114 4.1923 4.19 4.1894 4.1785 4.178 4.1694 4.1694 4.1694 4.1556 4.1498 4.1498 4.1357 4.1222 4.1222 4.1217 4.1192 4.1178 4.1139 4.1135 4.1125 4.1035 4.1025 4.1023 4.0971 4.0969 4.0915 4.0915 4.0914 4.0836 4.0804 4.0803 4.0722 4.065 4.065 4.0649 4.0644 4.0637 4.0616 4.0616 4.061 4.0572 4.0563 4.056 4.0545 4.0545 4.0522 4.0519 4.0514 4.0484 4.0467 4.0463 4.0422 4.0392 4.0388 4.0385 4.0385 4.0383 4.038 4.0379 4.0375 4.0364 4.0353 4.0344];
plot(1:100, curve)

---

编辑

我接受了Jonas给出的解决方案。基本上，对于曲线上的每个点p，我们找到与d 的最大距离为：

【问题讨论】：

我也想说画一条 45 度线：\

您希望您的图表与上述示例中的一般形状有多少偏差？换句话说，您是否期望图形的“肘部”总是靠近图形的同一个角？

@Amro：上半部分曲线的 100 点中只有 6 个点可以吗？试试plot(1:100,curve,'.')

@gnovice：我给出的例子是常见的情况，但这并不意味着规则没有例外 :) 就像我解释的那样，曲线代表模型可能性之间的权衡与其复杂性相比，因此您可以想象一个没有硬“弯头”的形状，而是看起来更像一条扁平线..

事实上，肯尼莫顿的答案是正确的。通过使用 AIC 等，您已经在某种程度上纠正了模型的复杂性。您应该使用其中一个标准并选择最低的，或者尝试在模型复杂度与 原始 拟合优度的曲线上找到肘部，而不是调整后的拟合优度 -复杂性。

【参考方案1】：

首先，快速回顾一下微积分：每个图形的一阶导数f' 表示函数f 的变化率。二阶导数f'' 表示f' 的变化率。如果f'' 很小，则意味着图形正在以适度的速度改变方向。但如果f'' 很大，则表示图形正在快速改变方向。

您想要隔离f'' 在图形域中最大的点。这些将是为您的最佳模型选择的候选点。您选择哪一点取决于您自己，因为您没有具体说明您对适合性与复杂性的重视程度。

【讨论】：

这个想法肯定是有效的，但最初的问题仍然存在；你如何指定阈值，之后你决定变化率太慢或太快。正如我之前描述的，我有大量的实验，这使得很难为所有情况设置一个通用值。

一般来说，你可以选择最大的f''。

这就是我们之前尝试的方式。然而，对于我们的应用程序来说，对有些嘈杂的数据进行两个导数证明不够稳健。

@Jonas 如果您有一个模型来适应模型复杂性的变化，这会更好。然后，您可以将噪声数据与模型拟合，并从模型中获取变化点。

@EsbenEickhardt 噪音当然会使任何评估变得复杂。您可以先尝试平滑。

【参考方案2】：

因此，解决此问题的一种方法是在肘部的 L 上放置两条线。但是由于曲线的一部分中只有几个点（正如我在评论中提到的那样），除非您检测到哪些点被间隔开并在它们之间插值以制造更统一的系列并且 ，否则线拟合会受到影响然后使用 RANSAC 找到适合 L 的两条线 - 有点复杂但并非不可能。

因此，这里有一个更简单的解决方案 - 由于 MATLAB 的缩放（显然），您创建的图表看起来像它们的方式。所以我所做的就是使用比例信息最小化从“原点”到您的点的距离。

请注意：原点估计可以显着改进，但我会留给你。

代码如下：

%% Order
curve = [8.4663 8.3457 5.4507 5.3275 4.8305 4.7895 4.6889 4.6833 4.6819 4.6542 4.6501 4.6287 4.6162 4.585 4.5535 4.5134 4.474 4.4089 4.3797 4.3494 4.3268 4.3218 4.3206 4.3206 4.3203 4.2975 4.2864 4.2821 4.2544 4.2288 4.2281 4.2265 4.2226 4.2206 4.2146 4.2144 4.2114 4.1923 4.19 4.1894 4.1785 4.178 4.1694 4.1694 4.1694 4.1556 4.1498 4.1498 4.1357 4.1222 4.1222 4.1217 4.1192 4.1178 4.1139 4.1135 4.1125 4.1035 4.1025 4.1023 4.0971 4.0969 4.0915 4.0915 4.0914 4.0836 4.0804 4.0803 4.0722 4.065 4.065 4.0649 4.0644 4.0637 4.0616 4.0616 4.061 4.0572 4.0563 4.056 4.0545 4.0545 4.0522 4.0519 4.0514 4.0484 4.0467 4.0463 4.0422 4.0392 4.0388 4.0385 4.0385 4.0383 4.038 4.0379 4.0375 4.0364 4.0353 4.0344];
x_axis = 1:numel(curve);
points = [x_axis ; curve ]'; %' - SO formatting

%% Get the scaling info
f = figure(1);
plot(points(:,1),points(:,2));
ticks = get(get(f,'CurrentAxes'),'YTickLabel');
ticks = str2num(ticks);
aspect = get(get(f,'CurrentAxes'),'DataAspectRatio');
aspect = [aspect(2) aspect(1)];    
close(f);   

%% Get the "origin"
O = [x_axis(1) ticks(1)];

%% Scale the data - now the scaled values look like MATLAB''s idea of
% what a good plot should look like
scaled_O = O.*aspect;
scaled_points = bsxfun(@times,points,aspect);

%% Find the closest point
del = sum((bsxfun(@minus,scaled_points,scaled_O).^2),2);
[val ind] = min(del);
best_ROC = [ind curve(ind)];

%% Display
plot(x_axis,curve,'.-');
hold on;
plot(O(1),O(2),'r*');
plot(best_ROC(1),best_ROC(2),'k*');

结果：

ALSO 对于Fit(maximize) 曲线，您必须将原点更改为[x_axis(1) ticks(end)]。

【讨论】：

一个有趣的想法，我唯一的问题是你必须绘制并使用 MATLAB 的自动缩放来找到原点。这不会很好地缩放，特别是因为我有超过一百万条曲线即时处理。 xaxis 总是保证从 1 开始，但我在 yaxis 上没有界限......并且回答你的问题，形状几乎是任意的，尽管我希望在曲线的开头快速上升/下降

所以我想问题是要弄清楚MATLAB的自动缩放......会检查出来。

【参考方案3】：

信息论模型选择的关键在于它已经考虑了参数的数量。因此，无需找到肘部，您只需找到最小值即可。

查找曲线的肘部仅在使用 fit 时相关。即使这样，您选择选择肘部的方法在某种意义上也是对参数数量的惩罚。要选择弯头，您需要最小化从原点到曲线的距离。距离计算中两个维度的相对权重会产生一个固有的惩罚项。信息论准则根据用于估计模型的参数数量和数据样本数量来设置此指标。

底线建议：使用 BIC 并取最小值。

【讨论】：

但是找到最小值并不是最优的。 IE。如果您查看 BIC 曲线，则会为位置 100 处的元素计算最小值。但是 20 和 100 的复杂度之间的差异非常小 - 模型的进一步复杂度的增益太小。

所以你是说 BIC 对模型复杂性的惩罚不足（我是真的，尽管我认为不是）。我的论点是，通过某种方法选择曲线的肘部会产生一个临时的未知惩罚项。如果您不喜欢 BIC 或 AIC 或任何其他基于 IC 的方法提供的答案，您最好制定一个惩罚项并使用它。只是我的意见。

在某种程度上我同意你的看法，因为在使用 AIC 和 BIC 时，我们没有比较绝对基线，而是以相对方式使用它们来比较模型（当所有其他事情是平等的）并且通过寻找这个 L 形，我们有效地引入了一个可能不平衡比赛场地的临时惩罚项。唯一的问题是选择 100 是因为计算成本和真正的上限因为复杂度接近 10000

这是真正的正确答案。通过使用 AIC 等，您已经在某种程度上纠正了模型的复杂性。您应该使用其中一个标准并选择最低的，或者尝试在模型复杂度与 原始 拟合优度而不是拟合优度的曲线上找到肘部 -调整复杂性。

【参考方案4】：

找到肘部的快速方法是从曲线的第一个点到最后一个点画一条线，然后找到离该线最远的数据点。

这当然在一定程度上取决于你在直线的平坦部分拥有的点数，但如果你每次测试相同数量的参数，它应该会相当不错。

curve = [8.4663 8.3457 5.4507 5.3275 4.8305 4.7895 4.6889 4.6833 4.6819 4.6542 4.6501 4.6287 4.6162 4.585 4.5535 4.5134 4.474 4.4089 4.3797 4.3494 4.3268 4.3218 4.3206 4.3206 4.3203 4.2975 4.2864 4.2821 4.2544 4.2288 4.2281 4.2265 4.2226 4.2206 4.2146 4.2144 4.2114 4.1923 4.19 4.1894 4.1785 4.178 4.1694 4.1694 4.1694 4.1556 4.1498 4.1498 4.1357 4.1222 4.1222 4.1217 4.1192 4.1178 4.1139 4.1135 4.1125 4.1035 4.1025 4.1023 4.0971 4.0969 4.0915 4.0915 4.0914 4.0836 4.0804 4.0803 4.0722 4.065 4.065 4.0649 4.0644 4.0637 4.0616 4.0616 4.061 4.0572 4.0563 4.056 4.0545 4.0545 4.0522 4.0519 4.0514 4.0484 4.0467 4.0463 4.0422 4.0392 4.0388 4.0385 4.0385 4.0383 4.038 4.0379 4.0375 4.0364 4.0353 4.0344];

%# get coordinates of all the points
nPoints = length(curve);
allCoord = [1:nPoints;curve]';              %'# SO formatting

%# pull out first point
firstPoint = allCoord(1,:);

%# get vector between first and last point - this is the line
lineVec = allCoord(end,:) - firstPoint;

%# normalize the line vector
lineVecN = lineVec / sqrt(sum(lineVec.^2));

%# find the distance from each point to the line:
%# vector between all points and first point
vecFromFirst = bsxfun(@minus, allCoord, firstPoint);

%# To calculate the distance to the line, we split vecFromFirst into two 
%# components, one that is parallel to the line and one that is perpendicular 
%# Then, we take the norm of the part that is perpendicular to the line and 
%# get the distance.
%# We find the vector parallel to the line by projecting vecFromFirst onto 
%# the line. The perpendicular vector is vecFromFirst - vecFromFirstParallel
%# We project vecFromFirst by taking the scalar product of the vector with 
%# the unit vector that points in the direction of the line (this gives us 
%# the length of the projection of vecFromFirst onto the line). If we 
%# multiply the scalar product by the unit vector, we have vecFromFirstParallel
scalarProduct = dot(vecFromFirst, repmat(lineVecN,nPoints,1), 2);
vecFromFirstParallel = scalarProduct * lineVecN;
vecToLine = vecFromFirst - vecFromFirstParallel;

%# distance to line is the norm of vecToLine
distToLine = sqrt(sum(vecToLine.^2,2));

%# plot the distance to the line
figure('Name','distance from curve to line'), plot(distToLine)

%# now all you need is to find the maximum
[maxDist,idxOfBestPoint] = max(distToLine);

%# plot
figure, plot(curve)
hold on
plot(allCoord(idxOfBestPoint,1), allCoord(idxOfBestPoint,2), 'or')

【讨论】：

谢谢我真的很喜欢这个解决方案！我不得不承认我很难理解你是如何计算点线距离的？

我同意这不是一个很好的评论。我试图更好地描述它。谁能想到几何最终会派上用场？

感谢您的解释。我最终发现这个页面刷新了我的几何：en.wikipedia.org/wiki/Vector_projection

如果您正在为相同的问题（我曾经）寻找 R 解决方案，它是 here。

这个conference paper, titled 'Finding a "Kneedle" in a Haystack: Detecting Knee Points in System Behavior'支持这个方法。

【参考方案5】：

我们可以用一种简单直观的方式说

如果我们从曲线上的任何一点到曲线的两个端点画两条线，这两条线以度为单位形成最小角度的点就是所需的点。

在这里，两条线可以可视化为手臂，点可视为肘点！

【参考方案6】：

如果有人需要上面 Jonas 发布的 Matlab 代码的工作 Python 版本。

import numpy as np
curve = [8.4663, 8.3457, 5.4507, 5.3275, 4.8305, 4.7895, 4.6889, 4.6833, 4.6819, 4.6542, 4.6501, 4.6287, 4.6162, 4.585, 4.5535, 4.5134, 4.474, 4.4089, 4.3797, 4.3494, 4.3268, 4.3218, 4.3206, 4.3206, 4.3203, 4.2975, 4.2864, 4.2821, 4.2544, 4.2288, 4.2281, 4.2265, 4.2226, 4.2206, 4.2146, 4.2144, 4.2114, 4.1923, 4.19, 4.1894, 4.1785, 4.178, 4.1694, 4.1694, 4.1694, 4.1556, 4.1498, 4.1498, 4.1357, 4.1222, 4.1222, 4.1217, 4.1192, 4.1178, 4.1139, 4.1135, 4.1125, 4.1035, 4.1025, 4.1023, 4.0971, 4.0969, 4.0915, 4.0915, 4.0914, 4.0836, 4.0804, 4.0803, 4.0722, 4.065, 4.065, 4.0649, 4.0644, 4.0637, 4.0616, 4.0616, 4.061, 4.0572, 4.0563, 4.056, 4.0545, 4.0545, 4.0522, 4.0519, 4.0514, 4.0484, 4.0467, 4.0463, 4.0422, 4.0392, 4.0388, 4.0385, 4.0385, 4.0383, 4.038, 4.0379, 4.0375, 4.0364, 4.0353, 4.0344]
nPoints = len(curve)
allCoord = np.vstack((range(nPoints), curve)).T
np.array([range(nPoints), curve])
firstPoint = allCoord[0]
lineVec = allCoord[-1] - allCoord[0]
lineVecNorm = lineVec / np.sqrt(np.sum(lineVec**2))
vecFromFirst = allCoord - firstPoint
scalarProduct = np.sum(vecFromFirst * np.matlib.repmat(lineVecNorm, nPoints, 1), axis=1)
vecFromFirstParallel = np.outer(scalarProduct, lineVecNorm)
vecToLine = vecFromFirst - vecFromFirstParallel
distToLine = np.sqrt(np.sum(vecToLine ** 2, axis=1))
idxOfBestPoint = np.argmax(distToLine)

【讨论】：

当我遍历“曲线”字典时，由于某种原因，这会失败，但在单次运行中表现完美。我不确定为什么循环会导致它失败，但运行时结果会有所不同，并且通常默认为“98”。

@user3486773 你检查过np.array([range(nPoints), curve])的代码行吗？

【参考方案7】：

双重派生方法。但是，它似乎不适用于嘈杂的数据。对于输出，您只需找到 d2 的最大值即可识别肘部。此实现在 R 中。

elbow_finder <- function(x_values, y_values) 
  i_max <- length(x_values) - 1

  # First and second derived
  first_derived <- list()
  second_derived <- list()

  # First derived
  for(i in 2:i_max)
    slope1 <- (y_values[i+1] - y_values[i]) / (x_values[i+1] - x_values[i])
    slope2 <- (y_values[i] - y_values[i-1]) / (x_values[i] - x_values[i-1])
    slope_avg <- (slope1 + slope2) / 2
    first_derived[[i]] <- slope_avg 
  
  first_derived[[1]] <- NA
  first_derived[[i_max+1]] <- NA
  first_derived <- unlist(first_derived)

  # Second derived
  for(i in 3:i_max-1)
    d1 <- (first_derived[i+1] - first_derived[i]) / (x_values[i+1] - x_values[i])
    d2 <- (first_derived[i] - first_derived[i-1]) / (x_values[i] - x_values[i-1])
    d_avg <- (d1 + d2) / 2
    second_derived[[i]] <- d_avg 
  
  second_derived[[1]] <- NA
  second_derived[[2]] <- NA
  second_derived[[i_max]] <- NA
  second_derived[[i_max+1]] <- NA
  second_derived <- unlist(second_derived)

  return(list(d1 = first_derived, d2 = second_derived))

【参考方案8】：

这是 Jonas 在 R 中实现的解决方案：

elbow_finder <- function(x_values, y_values) 
  # Max values to create line
  max_x_x <- max(x_values)
  max_x_y <- y_values[which.max(x_values)]
  max_y_y <- max(y_values)
  max_y_x <- x_values[which.max(y_values)]
  max_df <- data.frame(x = c(max_y_x, max_x_x), y = c(max_y_y, max_x_y))

  # Creating straight line between the max values
  fit <- lm(max_df$y ~ max_df$x)

  # Distance from point to line
  distances <- c()
  for(i in 1:length(x_values)) 
    distances <- c(distances, abs(coef(fit)[2]*x_values[i] - y_values[i] + coef(fit)[1]) / sqrt(coef(fit)[2]^2 + 1^2))
  

  # Max distance point
  x_max_dist <- x_values[which.max(distances)]
  y_max_dist <- y_values[which.max(distances)]

  return(c(x_max_dist, y_max_dist))

【参考方案9】：

我从事膝/肘点检测已经有一段时间了。无论如何，我是专家。 一些可能与此问题相关的方法。

DFDT 代表动态一阶导数阈值。它计算一阶导数并使用阈值算法来检测膝/肘点。 DSDT 类似但使用二阶导数，我的评估表明它们具有相似的性能。

S 方法是 L 方法的扩展。 L 方法将两条直线拟合到您的曲线，两条线之间的交点是膝/肘点。通过循环整体点、拟合线并评估 MSE（均方误差）来找到最佳拟合。 S 方法拟合 3 条直线，这提高了准确性，但也需要更多的计算。

我的所有代码都可以在 GitHub 上公开获得。此外，此article 可以帮助您找到有关该主题的更多信息。它只有四页长，所以应该很容易阅读。您可以使用代码，如果您想讨论任何方法，请随时这样做。

【参考方案10】：

如果您愿意，我已将其翻译为 R 作为自己的练习（请原谅我未优化的编码风格）。 *应用它在 k-means 上找到最佳集群数 - 工作得很好。

elbow.point = function(x)
elbow.curve = c(x)
nPoints = length(elbow.curve);
allCoord = cbind(c(1:nPoints),c(elbow.curve))
# pull out first point
firstPoint = allCoord[1,]
# get vector between first and last point - this is the line
lineVec = allCoord[nPoints,] - firstPoint;
# normalize the line vector
lineVecN = lineVec / sqrt(sum(lineVec^2));
# find the distance from each point to the line:
# vector between all points and first point
vecFromFirst = lapply(c(1:nPoints), function(x)
  allCoord[x,] - firstPoint
)
vecFromFirst = do.call(rbind, vecFromFirst)
rep.row<-function(x,n)
  matrix(rep(x,each=n),nrow=n)

scalarProduct = matrix(nrow = nPoints, ncol = 2)
scalarProduct[,1] = vecFromFirst[,1] * rep.row(lineVecN,nPoints)[,1]
scalarProduct[,2] = vecFromFirst[,2] * rep.row(lineVecN,nPoints)[,2]
scalarProduct = as.matrix(rowSums(scalarProduct))
vecFromFirstParallel = matrix(nrow = nPoints, ncol = 2)
vecFromFirstParallel[,1] = scalarProduct * lineVecN[1]
vecFromFirstParallel[,2] = scalarProduct * lineVecN[2]
vecToLine = lapply(c(1:nPoints), function(x)
  vecFromFirst[x,] - vecFromFirstParallel[x,]
)
vecToLine = do.call(rbind, vecToLine)
# distance to line is the norm of vecToLine
distToLine = as.matrix(sqrt(rowSums(vecToLine^2)))
##
which.max(distToLine)

函数的输入 x 应该是带有您的值的列表/向量

【参考方案11】：

不要忽视模型选择的 k 折交叉验证，它是 AIC/BIC 的绝佳替代品。还要考虑您正在建模的基本情况，您可以使用领域知识来帮助选择模型。