Jenks Natural Breaks en Python: ¿Cómo encontrar el número óptimo de descansos?

Encontré esta implementación de Python del algoritmo Jenks Natural Breaks y pude hacer que se ejecute en mi máquina con Windows 7. Es bastante rápido y encuentra los descansos en poco tiempo, considerando el tamaño de mis geodatos. Antes de usar este algoritmo de agrupamiento para mis datos, estaba usando sklearn.clustering.KMeans (aquí) el algoritmo. El problema que tuve con KMeans fue encontrar el parámetro óptimo del valor de K, pero lo "resolví" lanzando el algoritmo para diferentes valores de K y usando sklearn.metrics.silhouette_score (aquí) para encontrar el mejor K.

Mi pregunta es: si le digo al algoritmo de Natural Breaks que busque 5 clases (esa sería la K), ¿cómo puedo estar seguro de que este es el número de clases que mejor coinciden con mis datos? ¿Cómo validar que estoy eligiendo el mejor número de descansos?

¡Gracias!

Jenks Natural Breaks funciona optimizando la bondad de ajuste de varianza, un valor de 0 a 1 donde 0 = sin ajuste y 1 = ajuste perfecto. La clave para seleccionar el número de clases es encontrar un equilibrio entre detectar diferencias y sobreajustar sus datos. Para determinar el número óptimo de clases, le sugiero que use un valor umbral de GVF que desee y use primero el número de clases que satisfaga este valor.

A continuación se muestra una función para calcular el ajuste de bondad de la variación dada una matriz de valores para 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 .8, entonces podría incrementar el número de clases hasta que el GVF esté satisfecho:

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