Las pausas naturales de Jenks en Python: ¿Cómo encontrar el número óptimo de pausas?

Question

He encontrado esta implementación de Python de la Pausas naturales de Jenks algoritmo y pude hacerlo funcionar en mi máquina Windows 7. Es bastante rápido y encuentra las rupturas en poco tiempo, teniendo en cuenta el tamaño de mis geodatos. Antes de usar este algoritmo de clustering para mis datos, estaba usando sklearn.clustering.KMeans (aquí) algoritmo. El problema que tuve con KMeans, fue encontrar el parámetro óptimo del valor K, pero lo "resolví" lanzando el algoritmo para diferentes valores de K y utilizando sklearn.metrics.silhouette_score (aquí) para encontrar el mejor K.

Mi pregunta es: si le digo al algoritmo de Rupturas Naturales que encuentre 5 clases (que serían las K), ¿cómo puedo estar seguro de que éste es el número de clases que mejor se ajusta a mis datos? ¿Cómo validar que estoy eligiendo el mejor número de rupturas?

Gracias.

Answer 1

Jenks Natural Breaks funciona optimizando la bondad del ajuste de la varianza, un valor de 0 a 1 donde 0 = No hay ajuste y 1 = Ajuste perfecto. La clave en la selección del número de clases es encontrar un equilibrio entre la detección de diferencias y el sobreajuste de los datos. Para determinar el número óptimo de clases, le sugiero que utilice un valor umbral de GVF que desee y utilice el número de clases que satisfaga este valor primero.

A continuación se muestra una función para calcular la bondad del ajuste de la varianza dada una matriz de valores a clasificar y el número de clases seleccionadas:

from jenks import jenks
import numpy as np
def goodness_of_variance_fit(array, classes):
    # get the break points
    classes = jenks(array, classes)

    # do the actual classification
    classified = np.array([classify(i, classes) for i in array])

    # max value of zones
    maxz = max(classified)

    # nested list of zone indices
    zone_indices = [[idx for idx, val in enumerate(classified) if zone + 1 == val] for zone in range(maxz)]

    # sum of squared deviations from array mean
    sdam = np.sum((array - array.mean()) ** 2)

    # sorted polygon stats
    array_sort = [np.array([array[index] for index in zone]) for zone in zone_indices]

    # sum of squared deviations of class means
    sdcm = sum([np.sum((classified - classified.mean()) ** 2) for classified in array_sort])

    # goodness of variance fit
    gvf = (sdam - sdcm) / sdam

    return gvf

def classify(value, breaks):
    for i in range(1, len(breaks)):
        if value < breaks[i]:
            return i
    return len(breaks) - 1

Por ejemplo, considere que decide que el GVF debe ser al menos de 0,8, entonces podría incrementar el número de clases hasta que se satisfaga el GVF:

gvf = 0.0
nclasses = 2
while gvf < .8:
    gvf = goodness_of_variance_fit(array, nclasses)
    nclasses += 1