Jenks Natural Breaks en Python: Comment trouver le nombre optimal de pauses?

J'ai trouvé cette implémentation Python de l' algorithme Jenks Natural Breaks et je pouvais la faire fonctionner sur ma machine Windows 7. Il est assez rapide et trouve les ruptures en peu de temps, compte tenu de la taille de mes géodonnées. Avant d'utiliser cet algorithme de clustering pour mes données, j'utilisais sklearn.clustering.KMeans (ici) l' algorithme. Le problème que j'ai eu avec KMeans était de trouver le paramètre de valeur K optimal, mais je l'ai "résolu" en lançant l'algorithme pour différentes valeurs K et en utilisant sklearn.metrics.silhouette_score (ici) pour trouver le meilleur K.

Ma question est: si je dis à l'algorithme Natural Breaks de trouver 5 classes (ce serait le K), comment puis-je être sûr que c'est le nombre de classes qui correspondent le mieux à mes données? Comment valider que je choisis le meilleur nombre de pauses?

Merci!

Jenks Natural Breaks fonctionne en optimisant la qualité de l'ajustement de la variance, une valeur de 0 à 1 où 0 = pas d'ajustement et 1 = ajustement parfait. La clé pour sélectionner le nombre de classes est de trouver un équilibre entre la détection des différences et le surajustement de vos données. Pour déterminer le nombre optimal de classes, je vous suggère d'utiliser une valeur GVF seuil que vous désirez et d'utiliser le nombre de classes qui satisfait cette valeur en premier.

Vous trouverez ci-dessous une fonction pour calculer la qualité de l'ajustement de la variance en fonction d'un tableau de valeurs à classer et du nombre de classes sélectionnées:

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

Par exemple, considérez que vous décidez que le GVF doit être d'au moins 0,8, alors vous pouvez incrémenter le nombre de classes jusqu'à ce que le GVF soit satisfait:

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