algorithm - manualmente - punto de equilibrio en dinero

Encontrar el mejor punto de equilibrio en una curva (9)

Supongamos que tengo algunos datos, para los cuales quiero adaptar un modelo parametrizado. Mi objetivo es encontrar el mejor valor para este parámetro de modelo.

Estoy haciendo una selección de modelo usando un criterio de tipo AIC / BIC / MDL que premia modelos con un error bajo pero también penaliza a los modelos con alta complejidad (estamos buscando la explicación más simple pero más convincente para estos datos, por así decirlo, a la Occam navaja de afeitar ).

Siguiendo lo anterior, este es un ejemplo del tipo de cosas que obtengo para tres criterios diferentes (dos deben minimizarse y uno debe maximizarse):

Visualmente, puede ver fácilmente la forma del codo y puede elegir un valor para el parámetro en algún lugar de esa región. El problema es que estoy haciendo esto para una gran cantidad de experimentos y necesito una forma de encontrar este valor sin intervención.

Mi primera intuición fue tratar de dibujar una línea a 45 grados desde la esquina y seguir moviéndolo hasta que se cruza con la curva, pero eso es más fácil decirlo que hacerlo :) También puede pasar por alto la región de interés si la curva está algo sesgada.

¿Alguna idea sobre cómo implementar esto o mejores ideas?

Aquí están las muestras necesarias para reproducir una de las parcelas anteriores:

EDITAR

Acepté la solución dada por Jonas . Básicamente, para cada punto p en la curva, encontramos el que tiene la distancia máxima d dada por:

Aquí está la solución dada por Jonas implementada en 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)) }

De una manera simple e intuitiva podemos decir que

Si dibujamos dos líneas desde cualquier punto de la curva hasta los dos puntos finales de la curva, el punto en el que estas dos líneas forman el ángulo más pequeño en grados es el punto deseado.

¡Aquí, las dos líneas se pueden visualizar como los brazos y el punto como el punto del codo!

El método derivado doble Sin embargo, no parece funcionar bien para datos ruidosos. Para la salida, simplemente encontrará el valor máximo de d2 para identificar el codo. Esta implementación está en 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)) }

El punto de la selección del modelo teórico de la información es que ya cuenta para la cantidad de parámetros. Por lo tanto, no es necesario encontrar un codo, solo necesita encontrar el mínimo.

Encontrar el codo de la curva solo es relevante cuando se usa el ajuste. Incluso entonces, el método que elija para seleccionar el codo en cierto sentido establece una penalización para la cantidad de parámetros. Para seleccionar el codo, querrás minimizar la distancia desde el origen hasta la curva. La ponderación relativa de las dos dimensiones en el cálculo de distancia creará un término de penalización inherente. El criterio de la teoría de la información establece esta métrica en función del número de parámetros y la cantidad de muestras de datos utilizadas para estimar el modelo.

La recomendación de la línea de fondo: use BIC y tome el mínimo.

En caso de que alguien necesite una versión en Python que funcione del código de Matlab publicado por arriba.

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)

Entonces, una forma de resolver esto sería dos líneas ajustadas dos líneas a la L de tu codo. Pero dado que solo hay unos pocos puntos en una porción de la curva (como mencioné en el comentario), el ajuste de línea da un golpe a menos que usted detecte qué puntos están espaciados e interpola entre ellos para fabricar una serie más uniforme y luego usa RANSAC para encontrar dos líneas que se ajusten a la L , un poco intrincadas pero no imposibles.

Así que aquí hay una solución más simple: los gráficos que has puesto se ven como lo hacen gracias a la escala de MATLAB (obviamente). Entonces, todo lo que hice fue minimizar la distancia desde el "origen" a sus puntos usando la información de la escala.

Tenga en cuenta: la estimación de origen se puede mejorar de manera espectacular, pero se lo dejo a usted.

Aquí está el código:

%% 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*'');

Resultados:

TAMBIÉN para la curva Fit(maximize) tendrá que cambiar a origen en [x_axis(1) ticks(end)] .

He estado trabajando en la detección de puntos de rodilla / codo durante un tiempo. De ninguna manera, soy un experto. Algunos métodos que pueden ser relevantes para este problema.

DFDT significa Umbral Dinámico de Primera Derivación. Calcula la primera derivada y utiliza un algoritmo de Umbral para detectar el punto de rodilla / codo. DSDT es similar pero usa la segunda derivada, mi evaluación muestra que tienen desempeños similares.

El método S es una extensión del método L. El método L ajusta dos líneas rectas a su curva, la intercepción entre las dos líneas es el punto de la rodilla / codo. El mejor ajuste se encuentra haciendo un bucle de puntos generales, ajustando las líneas y evaluando el MSE (Mean Square Error). El método S se adapta a 3 líneas rectas, esto mejora la precisión pero también requiere un poco más de cálculo.

Todo mi código está disponible públicamente en GitHub . Además, este article puede ayudarlo a encontrar más información sobre el tema. Solo tiene cuatro páginas, por lo que debería ser fácil de leer. Puede usar el código y, si desea analizar alguno de los métodos, siéntase libre de hacerlo.

Primero, una revisión rápida del cálculo: la primera derivada f'' de cada gráfico representa la velocidad a la que cambia la función f está graficando. La segunda derivada f'''' representa la velocidad a la que f'' está cambiando. Si f'''' es pequeño, significa que el gráfico está cambiando de dirección a un ritmo modesto. Pero si f'''' es grande, significa que el gráfico está cambiando rápidamente de dirección.

Desea aislar los puntos en los que f'''' es más grande sobre el dominio del gráfico. Estos serán puntos candidatos para seleccionar para su modelo óptimo. El punto que elijas tendrá que depender de ti, ya que no has especificado exactamente cuánto valoras la aptitud frente a la complejidad.

Una manera rápida de encontrar el codo es dibujar una línea desde el primer hasta el último punto de la curva y luego encontrar el punto de datos que está más alejado de esa línea.

Por supuesto, esto depende de la cantidad de puntos que tenga en la parte plana de la línea, pero si prueba el mismo número de parámetros cada vez, debería salir razonablemente bien.

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