Нахождение лучшей точки компромисса на кривой
Скажем, у меня есть некоторые данные, для которых я хочу разместить над ними параметризованную модель. Моя цель - найти лучшее значение для этого параметра модели.
Я делаю выбор модели, используя критерий типа AIC/ BIC/ MDL, который вознаграждает модели с низкой ошибкой, но также оштрафовывает модели с высокой сложностью (мы ищем самое простое, но наиболее убедительное объяснение этих данных, так сказать, а-ля Оккама бритва).
Следуя вышесказанному, это пример того, что я получаю по трем различным критериям (два должны быть сведены к минимуму, а один - к максимальному):
Визуально вы можете легко увидеть форму локтя и выбрать значение для параметра где-нибудь в этой области. Проблема в том, что я делаю это для большого количества экспериментов, и мне нужен способ найти это значение без вмешательства.
Моя первая интуиция состояла в том, чтобы попытаться нарисовать линию под углом 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)
РЕДАКТИРОВАТЬ
Я принял решение, данное Джонасом. В основном, для каждой точки p на кривой мы находим тот с максимальным расстоянием d дано:
10 ответов
Быстрый способ найти колено - провести линию от первой до последней точки кривой, а затем найти точку данных, которая находится дальше всего от этой линии.
Это, конечно, в некоторой степени зависит от количества точек, которые у вас есть в плоской части линии, но если вы каждый раз тестируете одинаковое количество параметров, все должно получиться достаточно хорошо.
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')
В случае, если кому-то нужна рабочая версия кода Matlab на Python, опубликованная Jonas выше.
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)
Вот решение, данное Джонасом, реализованное в 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))
}
Во-первых, быстрый обзор исчисления: первая производная f' каждого графика представляет скорость, с которой функция f быть в курсе событий меняется. Вторая производная f'' представляет скорость, с которой f' меняется. Если f'' маленький, это означает, что график меняет направление в умеренном темпе. Но если f'' большой, это означает, что график быстро меняет направление.
Вы хотите выделить точки, в которых f'' является наибольшим в области графа. Это будут баллы-кандидаты для выбора оптимальной модели. Какой момент вы выберете, будет зависеть от вас, поскольку вы точно не указали, насколько вы цените фитнес по сравнению со сложностью.
Точка выбора теоретико-информационной модели заключается в том, что она уже учитывает количество параметров. Поэтому не нужно искать локоть, нужно только найти минимум.
Нахождение колена кривой актуально только при использовании подгонки. Даже тогда метод, который вы выбираете, чтобы выбрать колено, в некотором смысле устанавливает штраф за количество параметров. Чтобы выбрать колено, вы хотите минимизировать расстояние от начала координат до кривой. Относительный вес двух измерений в расчете расстояния создаст неотъемлемый штрафной член. Информационно-теоретический критерий задает эту метрику на основе количества параметров и количества выборок данных, используемых для оценки модели.
Итоговая рекомендация: используйте BIC и возьмите минимум.
Так что одним из способов решения этой проблемы было бы две подходящие две линии к 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*');
Результаты:
ТАКЖЕ для Fit(maximize) кривая вам придется изменить в начало координат [x_axis(1) ticks(end)],
Двойной производный метод. Это, однако, кажется, не работает хорошо для шумных данных. Для вывода вы просто найдете максимальное значение 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))
}
В простой и интуитивно понятной форме мы можем сказать, что
Если мы проведем две линии из любой точки кривой в обе конечные точки кривой, точка, в которой эти две линии образуют наименьший угол в градусах, является желаемой точкой.
Здесь две линии можно визуализировать как руки, а точку - как точку локтя!
Я работал над определением точки колена / локтя в течение некоторого времени. Я ни в коем случае не эксперт. Некоторые методы, которые могут иметь отношение к этой проблеме.
DFDT обозначает Порог Динамического Первичного Производного. Он вычисляет первую производную и использует алгоритм Thresholding для определения точки колена / локтя. DSDT похож, но использует вторую производную, моя оценка показывает, что они имеют аналогичные характеристики.
S-метод является расширением L-метода. L-метод соответствует двум прямым линиям вашей кривой, перехват между двумя линиями является точкой колена / локтя. Наилучшее совпадение достигается путем зацикливания общих точек, выравнивания линий и оценки MSE (среднеквадратическая ошибка). S-метод подходит для 3 прямых линий, это повышает точность, но также требует дополнительных вычислений.
Весь мой код общедоступен на GitHub. Кроме того, эта статья может помочь вам найти больше информации по теме. Это всего четыре страницы, поэтому его должно быть легко читать. Вы можете использовать код, и если вы хотите обсудить любой из методов, не стесняйтесь делать это.
Если хотите, я перевел его на R как упражнение для себя (простите за мой неоптимизированный стиль кодирования). * Применял его, чтобы найти оптимальное количество кластеров на k-средних - работал довольно хорошо.
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 функции должен быть списком / вектором с вашими значениями
Не пренебрегайте перекрестной проверкой в k-кратном порядке для выбора модели, отличной альтернативы AIC/BIC. Также подумайте о базовой ситуации, которую вы моделируете, и вы можете использовать знания предметной области, чтобы помочь выбрать модель.