Information loss by histograms

Question

Recently I became curious about what I imagine to be an old problem: the fidelity of histograms to an underlying data set. CrossValidated has a number of questions on the subject of "optimal histograms" such as:

Calculating optimal number of bins in a histogram

But the answers looked more like heuristics than metrics for quantifying the information that is lost from summarizing data. Can anyone point me in the direction of theoretical treatments (or algorithms) of this problem? As an example of the the type of answer I am looking for, I submitted an answer below.

Answer 1

I had a chance to examine the faithful data set in R, to score some plausible metrics: a two-sample Kolmogorov-Smirnov test (distance and p-value), Kullback–Leibler divergence, and Mutual Information for a range of different numbers of bins.

An overview of the data:

The duration between Old Faithful geyser eruptions can be plotted as a histogram, but what is the trade off of #bins vs. fidelity to the data set?

A 4 bin histogram shows too much averaging

A 14 bin histograms shows the two main modes

A 34 bin histogram show more minor modes

And a 54 bin histogram shows all of the data at the measurement accuracy:

Results of different metrics:

Using the metrics of information loss we can plot how they change with an increasing number of bins:

Here, the KS metrics remain mostly unchanged above 27 bins, KL Divergence (in my implementation) is a useless metric, but mutual information best captures the increasing accuracy with increasing bin count.

Thoughts?

Teh codez:

library(entropy)
library(ggplot2)
library(dplyr)
library(reshape2)
data("faithful")
duration <- faithful$waiting
xrange <- seq(min(duration),max(duration),by=1)
h.all <- hist(duration,c(xrange[1]-0.5,xrange+0.5))
metrics <- function(nbins){
  h <- hist(duration,seq(min(duration)-0.5*(max(duration)-min(duration)+1)/nbins,max(duration)+0.5*(max(duration)-min(duration)+1)/nbins,length.out=1+nbins)) # must force bin breakpoints
  x <- unlist(mapply(FUN=function(m,c){rep(m,c)},h$mids,h$counts))
  # KS test?
  ks <- ks.test(x,duration,alternative='two.sided',exact=FALSE)
  # KL divergence?
  KLD <- KL.empirical(x,duration)
  # Mutual information?
  bin.ix.orig <- findInterval(duration,h.all$breaks,rightmost.closed = TRUE)
  bin.ix.new <- findInterval(duration,h$breaks,rightmost.closed = TRUE)
  y2d <- table(data.frame(bin.ix.new,bin.ix.orig))
  MI <- mi.empirical(y2d)
  return(data.frame(KS.distance=as.numeric(ks$statistic), 
                    KS.p.value= ks$p.value, 
                    KL.Divergence = KLD,
                    Mutual.Info = MI))
}
df.sweep <- data.frame(nbins=seq(3,60)) %>% group_by(nbins) %>% do(metrics(.$nbins)) %>% ungroup()
ggplot(melt(df.sweep,id.vars='nbins'),aes(x=nbins,y=value))+geom_point()+facet_wrap(~variable,scales='free_y')

Answer 2

The optimal method to generate a histogram is to ditch it entirely and use a KDE instead.

Histograms are an anachronism, and with the advent of easy software to generate kernel density estimators there is no longer any reason to ever use them. Thus, in the practical sense, the optimal method to generate a histogram is to ditch it entirely and use a KDE instead. With the latter method, the problem of bin selection is not required. There is a bandwidth parameter that is estimated from the data (e.g., via MLE) and so the method does not depend on any arbitrary choice by the user. For this reason, asking for the optimal selection of bins in a histogram is a bit like asking for the optimal horse-shoe for racing in the Daytona 500.

Answer 3

I gave a talk a few years ago looking at estimation of polychoric correlations (which few statisticians have heard about -- these are MLE estimates of correlations underlying two ordinal variables, assuming that these underlying variables come from a bivariate normal distribution, and that the observed ordinal variables are a result of coarsening the underlying normal variables). I was simulating some data from such bivairate distribution (slides 7--9) to demonstrate the point, with two two variables categorized in three and four categories, respectively, which is way fewer than what a histogram would typically give you (but is typical with social science data response scales). On slide 15 of, I report estimation results: the standard error of the correlation between the underlying normal variables was 0.0146, while the standard error of the polychoric correlations between the ordinal variables was 0.0207. They estimate the same quantity, so the difference/ratio of the standard errors is the degree of information loss, which is pretty much exactly 50%.