Нахождение лучшей точки компромисса на кривой

Question

Скажем, у меня есть некоторые данные, для которых я хочу разместить над ними параметризованную модель. Моя цель - найти лучшее значение для этого параметра модели.

Я делаю выбор модели, используя тип критерия 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)

---

EDIT

Я принял решение, данное Йоной . В основном, для каждой точки p на кривой мы находим точку с максимальным расстоянием d, определяемым как:

Answer 1

Быстрый способ найти колено - это нарисовать линию от первой до последней точки кривой, а затем найти точку данных, которая находится дальше всего от этой линии.

Это, конечно, в некоторой степени зависит от количества точек, которые у вас есть в плоской части линии, но если вы каждый раз тестируете одинаковое количество параметров, оно должно получиться достаточно хорошо.

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

Answer 2

Если кому-то нужна рабочая Python версия кода Matlab , отправленная 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)

Answer 3

Суть выбора теоретико-информационной модели заключается в том, что она уже учитывает количество параметров. Поэтому не нужно искать локоть, нужно только найти минимум.

Нахождение колена кривой актуально только при использовании подгонки. Даже тогда метод, который вы выбираете, чтобы выбрать колено, в некотором смысле устанавливает штраф за количество параметров. Чтобы выбрать колено, вы хотите минимизировать расстояние от начала координат до кривой. Относительный вес двух измерений в расчете расстояния создаст неотъемлемый штрафной член. Информационно-теоретический критерий устанавливает эту метрику на основе количества параметров и количества выборок данных, использованных для оценки модели.

Итоговая рекомендация: используйте BIC и возьмите минимум.

Answer 4

Во-первых, краткий обзор исчисления: первая производная f' каждого графика представляет скорость изменения функции f в графике. Вторая производная f'' представляет скорость, с которой f' изменяется. Если f'' мало, это означает, что график меняет направление в умеренном темпе. Но если f'' большое, это означает, что график быстро меняет направление.

Вы хотите выделить точки, в которых f'' является наибольшим в области графика. Это будут баллы-кандидаты для выбора оптимальной модели. Какой момент вы выберете, будет зависеть от вас, так как вы точно не указали, насколько вы цените фитнес по сравнению со сложностью.

Answer 5

Вот решение, данное Джонасом, реализованное в 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))
}

Answer 6

Так что одним из способов решения этой проблемы было бы две подгонки двух линий к 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)].

Answer 7

В простой и интуитивно понятной форме мы можем сказать, что

Если мы рисуем две линии от любой точки кривой до обеих конечных точек кривой, точка, в которой эти две линии образуют наименьший угол в градусах, является желаемой точкой.

Здесь две линии можно представить как руки, а точку - как точку локтя!

Answer 8

Я работал над определением точки колена / локтя в течение некоторого времени. Я ни в коем случае не эксперт. Некоторые методы, которые могут иметь отношение к этой проблеме.

DFDT обозначает Порог Динамического Первичного Производного. Он вычисляет первую производную и использует алгоритм Thresholding для определения точки колена / локтя. DSDT похож, но использует вторую производную, моя оценка показывает, что они имеют схожие характеристики.

S-метод является расширением L-метода. L-метод соответствует двум прямым линиям вашей кривой, перехват между двумя линиями является точкой колена / локтя. Наилучшее совпадение достигается путем зацикливания общих точек, выравнивания линий и оценки MSE (среднеквадратическая ошибка) S-метод подходит для 3 прямых линий, это повышает точность, но также требует дополнительных вычислений.

Весь мой код общедоступен на GitHub . Кроме того, эта статья может помочь вам найти больше информации по теме. Это всего четыре страницы, поэтому его должно быть легко читать. Вы можете использовать код, и если вы хотите обсудить любой из методов, не стесняйтесь делать это.

Answer 9

Двойной производный метод. Это, однако, кажется, не работает хорошо для шумных данных. Для вывода вы просто найдете максимальное значение 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))
}

Answer 10

Не пренебрегайте k-кратной перекрестной проверкой для выбора модели, отличной альтернативы AIC / BIC Также подумайте о базовой ситуации, которую вы моделируете, и вам позволено использовать знания предметной области, чтобы помочь выбрать модель.