Quantitative thematic maps are typically classified based on the statistical distribution of the dataset. So, the first thing to do is create a histogram of the dataset. Depending on the shape of the distribution in the histogram, one selects a classification method. For example, the following typical rules apply for different types of distributions:
Multimodal distibution -> natural breaks
Rectangular distribution -> equal interval or quantile.
Normal distribution -> standard deviation
J-shaped distributions -> geometric progression
It has been argued that these classification schemes are based on the properties of the data’s statistical distribution without regard to the data’s spatial properties, particularly, the spatial autocorrelation in the data. There has been research into combining the statistical properties of the dataset with the spatial properties of the dataset. Unfortunately, many of these more complex data classification techniques have not become popular.
Despite their popularity, I find thematic maps to be a highly unstable way of representing reality. The picture created by these maps can instantly change depending on the classification technique used or the number of classes selected.
Looking at the histogram of datasets, I am amazed at how frequently boundary errors occur. These errors occur when the boundaries between the classed areas on the map do not align with the major breaks in the dataset. Ideally, the classification process should result in the boundaries created where there are minor breaks in the surface. In ArcGIS, one can always tweak the boundaries if you are not entirely satisfied with the statistical results, but this should be reported.