Jenks Natural Breaks in Python: Jak znaleźć optymalną liczbę przerw?

Znalazłem tej implementacji Pythona w Breaks Jenks Natural algorytmu i mogłem zrobić to uruchomić na moim Windows 7 maszynie. Jest dość szybki i wykrywa przerwy w krótkim czasie, biorąc pod uwagę rozmiar moich geodanych. Przed użyciem tego algorytmu grupowania dla moich danych korzystałem sklearn.clustering.KMeans (tutaj) z algorytmu. Problemem z KMeans było znalezienie optymalnego parametru wartości K, ale „rozwiązałem” to uruchomienie algorytmu dla różnych wartości K i użycie sklearn.metrics.silhouette_score (tutaj) do znalezienia najlepszego K.

Moje pytanie brzmi: jeśli powiem algorytmowi Natural Breaks, aby znalazł 5 klas (to byłby K), jak mogę się upewnić, że jest to liczba klas, które najlepiej pasują do moich danych? Jak potwierdzić, że wybieram najlepszą liczbę przerw?

Dzięki!

Jenks Natural Breaks działa poprzez optymalizację dobroci wariancji dopasowania, wartość od 0 do 1, gdzie 0 = brak dopasowania i 1 = idealne dopasowanie. Kluczem do wyboru liczby klas jest znalezienie równowagi między wykrywaniem różnic a nadmiernym przetwarzaniem danych. Aby określić optymalną liczbę klas, sugeruję użycie pożądanej wartości progowej GVF i użycie liczby klas, które spełniają tę wartość w pierwszej kolejności.

Poniżej znajduje się funkcja do obliczania dobroci dopasowania wariancji, biorąc pod uwagę tablicę wartości do sklasyfikowania i liczbę wybranych klas:

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

Na przykład, rozważ, że zdecydujesz, że GVF powinien wynosić co najmniej 0,8, wtedy możesz zwiększać liczbę klas, aż GVF będzie spełniony:

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