Jenks Natural Breaks in Python: How to find the optimum number of breaks?

I found this Python implementation of the Jenks Natural Breaks algorithm and I could make it run on my Windows 7 machine. It is pretty fast and it finds the breaks in few time, considering the size of my geodata. Before using this clustering algorithm for my data, I was using sklearn.clustering.KMeans (here) algorithm. The problem I had with KMeans, was finding the optimum K value parameter, but I "solved" it launching the algorithm for different K values and using sklearn.metrics.silhouette_score (here) to find the best K.

My question is: if I tell the Natural Breaks algorithm to find 5 classes (that would be the K), how can I be sure that this is the number of classes that best match my data? How to validate that I am choosing the best number of breaks?

Thanks!

• So that we can objectively determine what "best" means, could you explain the sense in which the classes "match" data? (Or, really, how you would quantify any degree of mis-matching.)
– whuber
Commented Mar 29, 2015 at 20:45
• Using Silhouette with Jenks should be comparable to using it with kmeans. It's a heuristic and you shouldn't trust it blindly. IMHO the best is to visualize your rrsults. Commented Mar 30, 2015 at 5:56
• Whuber: Best, using Silhouette, would mean the number of classes that make the index be closer to 1, according to the definition on sklearn site: scikit-learn.org/stable/modules/generated/… Anony-Mousse: I cannot visualize 20+ variables, prepare maps for that and expect my brain not to mess up with the number of classes. I need to rely on an index that says, "for variable X, the best you can do is to use Y classes". Moreover I need to re-run the analysis several times, the viz approach is slow unfortunately... Commented Mar 30, 2015 at 8:02
• from jenks import jenks: gives following error Traceback (most recent call last): File "<stdin>", line 1, in <module> ImportError: cannot import name jenks Commented Jun 23, 2016 at 6:08

Jenks Natural Breaks works by optimizing the Goodness of Variance Fit, a value from 0 to 1 where 0 = No Fit and 1 = Perfect Fit. The key in selecting the number of classes is to find a balance between detecting differences and overfitting your data. To determine the optimum number of classes, I suggest you use a threshold GVF value you desire and use the number of classes that satisfies this value first.

Below is a function to calculate the Goodness of Variance Fit given an array of values to classify and the number of classes selected:

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


For example, consider you decide the GVF should be at least .8, then you could increment the number of classes until the GVF is satisfied:

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

• Hi, useful post. Can I know how to interpret Gvf = 87% for n_classes = 4? Does it mean with 4 classes, we are able to seperate the data (into bins) correctly by 87%? Commented Sep 25, 2022 at 3:38