diff --git a/examples/statistics/boxplot_color_demo.py b/examples/statistics/boxplot_color_demo.py
index 474f427d76eb..085a5a9e64c5 100644
--- a/examples/statistics/boxplot_color_demo.py
+++ b/examples/statistics/boxplot_color_demo.py
@@ -1,10 +1,16 @@
-"""Box plots with custom fill colors.
+"""
+=================================
+Box plots with custom fill colors
+=================================
 
 This plot illustrates how to create two types of box plots
 (rectangular and notched), and how to fill them with custom
 colors by accessing the properties of the artists of the
 box plots. Additionally, the ``labels`` parameter is used to
 provide x-tick labels for each sample.
+
+A good general reference on boxplots and their history can be found
+here: http://vita.had.co.nz/papers/boxplots.pdf
 """
 
 import matplotlib.pyplot as plt
@@ -15,7 +21,7 @@
 all_data = [np.random.normal(0, std, size=100) for std in range(1, 4)]
 labels = ['x1', 'x2', 'x3']
 
-fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
+fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4))
 
 # rectangular box plot
 bplot1 = axes[0].boxplot(all_data,
diff --git a/examples/statistics/boxplot_demo.py b/examples/statistics/boxplot_demo.py
index 19849469b83d..d065fd5997cb 100644
--- a/examples/statistics/boxplot_demo.py
+++ b/examples/statistics/boxplot_demo.py
@@ -1,5 +1,19 @@
 """
-Demo of the new boxplot functionality
+=========================================
+Demo of artist customization in box plots
+=========================================
+
+This example demonstrates how to use the various kwargs
+to fully customize box plots. The first figure demonstrates
+how to remove and add individual components (note that the
+mean is the only value not shown by default). The second
+figure demonstrates how the styles of the artists can
+be customized. It also demonstrates how to set the limit
+of the whiskers to specific percentiles (lower right axes)
+
+A good general reference on boxplots and their history can be found
+here: http://vita.had.co.nz/papers/boxplots.pdf
+
 """
 
 import numpy as np
@@ -23,7 +37,8 @@
 axes[0, 2].set_title('showmeans=True,\nmeanline=True', fontsize=fs)
 
 axes[1, 0].boxplot(data, labels=labels, showbox=False, showcaps=False)
-axes[1, 0].set_title('Tufte Style \n(showbox=False,\nshowcaps=False)', fontsize=fs)
+tufte_title = 'Tufte Style \n(showbox=False,\nshowcaps=False)'
+axes[1, 0].set_title(tufte_title, fontsize=fs)
 
 axes[1, 1].boxplot(data, labels=labels, notch=True, bootstrap=10000)
 axes[1, 1].set_title('notch=True,\nbootstrap=10000', fontsize=fs)
@@ -62,7 +77,8 @@
                    showmeans=True)
 axes[1, 0].set_title('Custom mean\nas point', fontsize=fs)
 
-axes[1, 1].boxplot(data, meanprops=meanlineprops, meanline=True, showmeans=True)
+axes[1, 1].boxplot(data, meanprops=meanlineprops, meanline=True,
+                   showmeans=True)
 axes[1, 1].set_title('Custom mean\nas line', fontsize=fs)
 
 axes[1, 2].boxplot(data, whis=[15, 85])
diff --git a/examples/statistics/boxplot_vs_violin_demo.py b/examples/statistics/boxplot_vs_violin_demo.py
index c4150e5db7bb..f183a6c7487d 100644
--- a/examples/statistics/boxplot_vs_violin_demo.py
+++ b/examples/statistics/boxplot_vs_violin_demo.py
@@ -1,19 +1,29 @@
-"""Box plot - violin plot comparison.
+"""
+===================================
+Box plot vs. violin plot comparison
+===================================
+
+Note that although violin plots are closely related to Tukey's (1977)
+box plots, they add useful information such as the distribution of the
+sample data (density trace).
 
-Note that although violin plots are closely related to Tukey's (1977) box
-plots, they add useful information such as the distribution of the sample
-data (density trace).
+By default, box plots show data points outside 1.5 * the inter-quartile
+range as outliers above or below the whiskers whereas violin plots show
+the whole range of the data.
 
-By default, box plots show data points outside 1.5 x the inter-quartile range
-as outliers above or below the whiskers whereas violin plots show the whole
-range of the data.
+A good general reference on boxplots and their history can be found
+here: http://vita.had.co.nz/papers/boxplots.pdf
 
 Violin plots require matplotlib >= 1.4.
+
+For more information on violin plots, the scikit-learn docs have a great
+section: http://scikit-learn.org/stable/modules/density.html
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
-fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
+fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4))
 
 # generate some random test data
 all_data = [np.random.normal(0, std, 100) for std in range(6, 10)]
diff --git a/examples/statistics/bxp_demo.py b/examples/statistics/bxp_demo.py
index 85d21775eca8..90d8b4c4e399 100644
--- a/examples/statistics/bxp_demo.py
+++ b/examples/statistics/bxp_demo.py
@@ -1,5 +1,17 @@
 """
-Demo of the new boxplot drawer function
+===================================
+Demo of the boxplot drawer function
+===================================
+
+This example demonstrates how to pass pre-computed box plot
+statistics to the box plot drawer. The first figure demonstrates
+how to remove and add individual components (note that the
+mean is the only value not shown by default). The second
+figure demonstrates how the styles of the artists can
+be customized.
+
+A good general reference on boxplots and their history can be found
+here: http://vita.had.co.nz/papers/boxplots.pdf
 """
 
 import numpy as np
@@ -21,6 +33,7 @@
     stats[n]['mean'] *= 2
 
 print(stats[0].keys())
+
 fs = 10  # fontsize
 
 # demonstrate how to toggle the display of different elements:
@@ -35,7 +48,8 @@
 axes[0, 2].set_title('showmeans=True,\nmeanline=True', fontsize=fs)
 
 axes[1, 0].bxp(stats, showbox=False, showcaps=False)
-axes[1, 0].set_title('Tufte Style\n(showbox=False,\nshowcaps=False)', fontsize=fs)
+tufte_title = 'Tufte Style\n(showbox=False,\nshowcaps=False)'
+axes[1, 0].set_title(tufte_title, fontsize=fs)
 
 axes[1, 1].bxp(stats, shownotches=True)
 axes[1, 1].set_title('notch=True', fontsize=fs)
@@ -50,7 +64,6 @@
 fig.subplots_adjust(hspace=0.4)
 plt.show()
 
-
 # demonstrate how to customize the display different elements:
 boxprops = dict(linestyle='--', linewidth=3, color='darkgoldenrod')
 flierprops = dict(marker='o', markerfacecolor='green', markersize=12,
@@ -71,7 +84,8 @@
                showmeans=True)
 axes[1, 0].set_title('Custom mean\nas point', fontsize=fs)
 
-axes[1, 1].bxp(stats, meanprops=meanlineprops, meanline=True, showmeans=True)
+axes[1, 1].bxp(stats, meanprops=meanlineprops, meanline=True,
+               showmeans=True)
 axes[1, 1].set_title('Custom mean\nas line', fontsize=fs)
 
 for ax in axes.flatten():
diff --git a/examples/statistics/customized_violin_demo.py b/examples/statistics/customized_violin_demo.py
index 260842a2fed8..29dfda7b894e 100644
--- a/examples/statistics/customized_violin_demo.py
+++ b/examples/statistics/customized_violin_demo.py
@@ -1,6 +1,18 @@
-# Customizing violin plots
-#
-#
+"""
+=================================
+Demo of violin plot customization
+=================================
+
+This example demonstrates how to fully customize violin plots.
+The first plot shows the default style by providing only
+the data. The second plot first limits what matplotlib draws
+with additional kwargs. Then a simplified representation of
+a box plot is drawn on top. Lastly, the styles of the artists
+of the violins are modified.
+
+For more information on violin plots, the scikit-learn docs have a great
+section: http://scikit-learn.org/stable/modules/density.html
+"""
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -22,21 +34,28 @@ def percentile(vals, p):
 def adjacent_values(vals):
     q1 = percentile(vals, 0.25)
     q3 = percentile(vals, 0.75)
-    uav = q3 + (q3-q1)*1.5
+    iqr = q3 - q1  # inter-quartile range
+
+    # upper adjacent values
+    uav = q3 + iqr * 1.5
     if uav > vals[-1]:
         uav = vals[-1]
     if uav < q3:
         uav = q3
-    lav = q1 - (q3-q1)*1.5
+
+    # lower adjacent values
+    lav = q1 - iqr * 1.5
     if lav < vals[0]:
         lav = vals[0]
     if lav > q1:
         lav = q1
     return [lav, uav]
 
+
 # create test data
-dat = [np.random.normal(0, std, 100) for std in range(6, 10)]
-lab = ['a', 'b', 'c', 'd']    # labels
+np.random.seed(123)
+dat = [np.random.normal(0, std, 100) for std in range(1, 5)]
+lab = ['A', 'B', 'C', 'D']    # labels
 med = []    # medians
 iqr = []    # inter-quantile ranges
 avs = []    # upper and lower adjacent values
@@ -47,15 +66,23 @@ def adjacent_values(vals):
     avs.append(adjacent_values(sarr))
 
 # plot the violins
-fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 5))
-parts = ax.violinplot(dat, showmeans=False, showmedians=False,
-                      showextrema=False)
+fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(9, 4),
+                               sharey=True)
+_ = ax1.violinplot(dat)
+parts = ax2.violinplot(dat, showmeans=False, showmedians=False,
+                       showextrema=False)
+
+ax1.set_title('Default violin plot')
+ax2.set_title('Customized violin plot')
 
-# plot medians and averages
+# plot whiskers as thin lines, quartiles as fat lines,
+# and medians as points
 for i in range(len(med)):
-    ax.plot([i+1, i+1], avs[i], '-', c='black', lw=1)
-    ax.plot([i+1, i+1], iqr[i], '-', c='black', lw=5)
-    ax.plot(i+1, med[i], 'o', mec='none', c='white', ms=6)
+    # whiskers
+    ax2.plot([i + 1, i + 1], avs[i], '-', color='black', linewidth=1)
+    ax2.plot([i + 1, i + 1], iqr[i], '-', color='black', linewidth=5)
+    ax2.plot(i + 1, med[i], 'o', color='white',
+             markersize=6, markeredgecolor='none')
 
 # customize colors
 for pc in parts['bodies']:
@@ -63,15 +90,15 @@ def adjacent_values(vals):
     pc.set_edgecolor('black')
     pc.set_alpha(1)
 
-ax.get_xaxis().set_tick_params(direction='out')
-ax.xaxis.set_ticks_position('bottom')
-ax.set_xticks([x+1 for x in range(len(lab))])
-ax.set_xticklabels(lab)
-ax.set_xlim(0.25, len(lab)+0.75)
-ax.set_ylabel('ylabel')
-ax.set_xlabel('xlabel')
-ax.set_title('customized violin plot')
+ax1.set_ylabel('Observed values')
+for ax in [ax1, ax2]:
+    ax.get_xaxis().set_tick_params(direction='out')
+    ax.xaxis.set_ticks_position('bottom')
+    ax.set_xticks(np.arange(1, len(lab) + 1))
+    ax.set_xticklabels(lab)
+    ax.set_xlim(0.25, len(lab) + 0.75)
+    ax.set_xlabel('Sample name')
 
-plt.subplots_adjust(bottom=0.15)
+plt.subplots_adjust(bottom=0.15, wspace=0.05)
 
 plt.show()
diff --git a/examples/statistics/errorbar_demo.py b/examples/statistics/errorbar_demo.py
index d970e62d8441..540cec3ab9b7 100644
--- a/examples/statistics/errorbar_demo.py
+++ b/examples/statistics/errorbar_demo.py
@@ -1,6 +1,13 @@
 """
-Demo of the errorbar function.
+=============================
+Demo of the errorbar function
+=============================
+
+This exhibits the most basic use of the error bar method.
+In this case, constant values are provided for the error
+in both the x- and y-directions.
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -8,5 +15,6 @@
 x = np.arange(0.1, 4, 0.5)
 y = np.exp(-x)
 
-plt.errorbar(x, y, xerr=0.2, yerr=0.4)
+fig, ax = plt.subplots()
+ax.errorbar(x, y, xerr=0.2, yerr=0.4)
 plt.show()
diff --git a/examples/statistics/errorbar_demo_features.py b/examples/statistics/errorbar_demo_features.py
index 4ec7a68ed437..dfe8f68b5524 100644
--- a/examples/statistics/errorbar_demo_features.py
+++ b/examples/statistics/errorbar_demo_features.py
@@ -1,37 +1,46 @@
 """
-Demo of errorbar function with different ways of specifying error bars.
+===================================================
+Demo of the different ways of specifying error bars
+===================================================
 
-Errors can be specified as a constant value (as shown in `errorbar_demo.py`),
-or as demonstrated in this example, they can be specified by an N x 1 or 2 x N,
-where N is the number of data points.
+Errors can be specified as a constant value (as shown in
+`errorbar_demo.py`). However, this example demonstrates
+how they vary by specifying arrays of error values.
 
-N x 1:
-    Error varies for each point, but the error values are symmetric (i.e. the
-    lower and upper values are equal).
+If the raw ``x`` and ``y`` data have length N, there are two options:
 
-2 x N:
-    Error varies for each point, and the lower and upper limits (in that order)
-    are different (asymmetric case)
+Array of shape (N,):
+    Error varies for each point, but the error values are
+    symmetric (i.e. the lower and upper values are equal).
 
-In addition, this example demonstrates how to use log scale with errorbar.
+Array of shape (2, N):
+    Error varies for each point, and the lower and upper limits
+    (in that order) are different (asymmetric case)
+
+In addition, this example demonstrates how to use log
+scale with error bars.
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
 # example data
 x = np.arange(0.1, 4, 0.5)
 y = np.exp(-x)
+
 # example error bar values that vary with x-position
 error = 0.1 + 0.2 * x
-# error bar values w/ different -/+ errors
-lower_error = 0.4 * error
-upper_error = error
-asymmetric_error = [lower_error, upper_error]
 
 fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True)
 ax0.errorbar(x, y, yerr=error, fmt='-o')
 ax0.set_title('variable, symmetric error')
 
+# error bar values w/ different -/+ errors that
+# also vary with the x-position
+lower_error = 0.4 * error
+upper_error = error
+asymmetric_error = [lower_error, upper_error]
+
 ax1.errorbar(x, y, xerr=asymmetric_error, fmt='o')
 ax1.set_title('variable, asymmetric error')
 ax1.set_yscale('log')
diff --git a/examples/statistics/errorbar_limits.py b/examples/statistics/errorbar_limits.py
index f12e0e3147d9..a765445598fb 100644
--- a/examples/statistics/errorbar_limits.py
+++ b/examples/statistics/errorbar_limits.py
@@ -1,50 +1,76 @@
 """
-Demo of the errorbar function, including upper and lower limits
+===========================================================
+Demo of how to include upper and lower limits in error bars
+===========================================================
+
+In matplotlib, errors bars can have "limits". Applying limits to the
+error bars essentially makes the error unidirectional. Because of that,
+upper and lower limits can be applied in both the y- and x-directions
+via the ``uplims``, ``lolims``, ``xuplims``, and ``xlolims`` parameters,
+respectively. These parameters can be scalar or boolean arrays.
+
+For example, if ``xlolims`` is ``True``, the x-error bars will only
+extend from the data towards increasing values. If ``uplims`` is an
+array filled with ``False`` except for the 4th and 7th values, all of the
+y-error bars will be bidirectional, except the 4th and 7th bars, which
+will extend from the data towards decreasing y-values.
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
 # example data
-x = np.arange(0.5, 5.5, 0.5)
+x = np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
 y = np.exp(-x)
 xerr = 0.1
 yerr = 0.2
+
+# lower & upper limits of the error
+lolims = np.array([0, 0, 1, 0, 1, 0, 0, 0, 1, 0], dtype=bool)
+uplims = np.array([0, 1, 0, 0, 0, 1, 0, 0, 0, 1], dtype=bool)
 ls = 'dotted'
 
-fig = plt.figure()
-ax = fig.add_subplot(1, 1, 1)
+fig, ax = plt.subplots(figsize=(7, 4))
 
 # standard error bars
-plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls=ls)
+ax.errorbar(x, y, xerr=xerr, yerr=yerr, linestyle=ls)
 
 # including upper limits
-uplims = np.zeros(x.shape)
-uplims[[1, 5, 9]] = True
-plt.errorbar(x, y + 0.5, xerr=xerr, yerr=yerr, uplims=uplims, ls=ls)
+ax.errorbar(x, y + 0.5, xerr=xerr, yerr=yerr, uplims=uplims,
+            linestyle=ls)
 
 # including lower limits
-lolims = np.zeros(x.shape)
-lolims[[2, 4, 8]] = True
-plt.errorbar(x, y + 1.0, xerr=xerr, yerr=yerr, lolims=lolims, ls=ls)
+ax.errorbar(x, y + 1.0, xerr=xerr, yerr=yerr, lolims=lolims,
+            linestyle=ls)
 
 # including upper and lower limits
-plt.errorbar(x, y + 1.5, marker='o', ms=8, xerr=xerr, yerr=yerr,
-             lolims=lolims, uplims=uplims, ls=ls)
+ax.errorbar(x, y + 1.5, xerr=xerr, yerr=yerr,
+            lolims=lolims, uplims=uplims,
+            marker='o', markersize=8,
+            linestyle=ls)
 
-# including xlower and xupper limits
+# Plot a series with lower and upper limits in both x & y
+# constant x-error with varying y-error
 xerr = 0.2
 yerr = np.zeros(x.shape) + 0.2
 yerr[[3, 6]] = 0.3
+
+# mock up some limits by modifying previous data
 xlolims = lolims
 xuplims = uplims
 lolims = np.zeros(x.shape)
 uplims = np.zeros(x.shape)
-lolims[[6]] = True
-uplims[[3]] = True
-plt.errorbar(x, y + 2.1, marker='o', ms=8, xerr=xerr, yerr=yerr,
-             xlolims=xlolims, xuplims=xuplims, uplims=uplims, lolims=lolims,
-             ls='none')
+lolims[[6]] = True  # only limited at this index
+uplims[[3]] = True  # only limited at this index
+
+# do the plotting
+ax.errorbar(x, y + 2.1, xerr=xerr, yerr=yerr,
+            xlolims=xlolims, xuplims=xuplims,
+            uplims=uplims, lolims=lolims,
+            marker='o', markersize=8,
+            linestyle='none')
 
+# tidy up the figure
 ax.set_xlim((0, 5.5))
 ax.set_title('Errorbar upper and lower limits')
 plt.show()
diff --git a/examples/statistics/errorbars_and_boxes.py b/examples/statistics/errorbars_and_boxes.py
index a0a421640c58..113e299be236 100644
--- a/examples/statistics/errorbars_and_boxes.py
+++ b/examples/statistics/errorbars_and_boxes.py
@@ -1,6 +1,25 @@
-'''
-Example to create boxes from error using PatchCollection
-'''
+"""
+============================================================
+Demo on creating boxes from error bars using PatchCollection
+============================================================
+
+In this example, we snazz up a pretty standard error bar plot by adding
+a rectangle patch defined by the limits of the bars in both the x- and
+y- directions. To do this, we have to write our own custom function
+called ``make_error_boxes``. Close inspection of this function will
+reveal the preferred pattern in writing functions for matplotlib:
+
+  1. an ``Axes`` object is passed directly to the function
+  2. the function operates on the `Axes` methods directly, not through
+     the ``pyplot`` interface
+  3. plotting kwargs that could be abbreviated are spelled out for
+     better code readability in the future (for example we use
+     ``facecolor`` instead of ``fc``)
+  4. the artists returned by the ``Axes`` plotting methods are then
+     returned by the function so that, if desired, their styles
+     can be modified later outside of the function (they are not
+     modified in this example).
+"""
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -11,38 +30,44 @@
 n = 5
 
 # Dummy data
+np.random.seed(10)
 x = np.arange(0, n, 1)
-y = np.random.rand(n)*5.
+y = np.random.rand(n) * 5.
 
 # Dummy errors (above and below)
-xerr = np.random.rand(2, n)
-yerr = np.random.rand(2, n)
-
-# Create figure and axes
-fig, ax = plt.subplots(1)
+xerr = np.random.rand(2, n) + 0.1
+yerr = np.random.rand(2, n) + 0.2
 
 
-def make_error_boxes(ax, xdata, ydata, xerror, yerror, fc='r', ec='None', alpha=0.5):
+def make_error_boxes(ax, xdata, ydata, xerror, yerror, facecolor='r',
+                     edgecolor='None', alpha=0.5):
 
     # Create list for all the error patches
     errorboxes = []
 
     # Loop over data points; create box from errors at each point
-    for xc, yc, xe, ye in zip(xdata, ydata, xerror.T, yerror.T):
-        rect = Rectangle((xc-xe[0], yc-ye[0]), xe.sum(), ye.sum())
+    for x, y, xe, ye in zip(xdata, ydata, xerror.T, yerror.T):
+        rect = Rectangle((x - xe[0], y - ye[0]), xe.sum(), ye.sum())
         errorboxes.append(rect)
 
     # Create patch collection with specified colour/alpha
-    pc = PatchCollection(errorboxes, facecolor=fc, alpha=alpha, edgecolor=ec)
+    pc = PatchCollection(errorboxes, facecolor=facecolor, alpha=alpha,
+                         edgecolor=edgecolor)
 
     # Add collection to axes
     ax.add_collection(pc)
 
     # Plot errorbars
-    ax.errorbar(xdata, ydata, xerr=xerror, yerr=yerror, fmt='None', ecolor='k')
+    artists = ax.errorbar(xdata, ydata, xerr=xerror, yerr=yerror,
+                          fmt='None', ecolor='k')
 
+    return artists
+
+
+# Create figure and axes
+fig, ax = plt.subplots(1)
 
 # Call function to create error boxes
-make_error_boxes(ax, x, y, xerr, yerr)
+_ = make_error_boxes(ax, x, y, xerr, yerr)
 
 plt.show()
diff --git a/examples/statistics/histogram_demo_cumulative.py b/examples/statistics/histogram_demo_cumulative.py
index 23b5c8c93463..e5918684e2f2 100644
--- a/examples/statistics/histogram_demo_cumulative.py
+++ b/examples/statistics/histogram_demo_cumulative.py
@@ -1,29 +1,72 @@
 """
-Demo of the histogram (hist) function used to plot a cumulative distribution.
+==========================================================
+Demo of using histograms to plot a cumulative distribution
+==========================================================
+
+This shows how to plot a cumulative, normalized histogram as a
+step function in order to visualize the empirical cumulative
+distribution function (CDF) of a sample. We also use the ``mlab``
+module to show the theoretical CDF.
+
+A couple of other options to the ``hist`` function are demonstrated.
+Namely, we use the ``normed`` parameter to normalize the histogram and
+a couple of different options to the ``cumulative`` parameter.
+The ``normed`` parameter takes a boolean value. When ``True``, the bin
+heights are scaled such that the total area of the histogram is 1. The
+``cumulative`` kwarg is a little more nuanced. Like ``normed``, you
+can pass it True or False, but you can also pass it -1 to reverse the
+distribution.
+
+Since we're showing a normalized and cumulative histogram, these curves
+are effectively the cumulative distribution functions (CDFs) of the
+samples. In engineering, empirical CDFs are sometimes called
+"non-exceedance" curves. In other words, you can look at the
+y-value for a given-x-value to get the probability of and observation
+from the sample not exceeding that x-value. For example, the value of
+225 on the x-axis corresponds to about 0.85 on the y-axis, so there's an
+85% chance that an observation in the sample does not exceed 225.
+Conversely, setting, ``cumulative`` to -1 as is done in the
+last series for this example, creates a "exceedance" curve.
+
+Selecting different bin counts and sizes can significantly affect the
+shape of a histogram. The Astropy docs have a great section on how to
+select these parameters:
+http://docs.astropy.org/en/stable/visualization/histogram.html
 
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib import mlab
 
+np.random.seed(0)
 
 mu = 200
 sigma = 25
 n_bins = 50
-x = mu + sigma*np.random.randn(10000)
+x = np.random.normal(mu, sigma, size=100)
+
+fig, ax = plt.subplots(figsize=(8, 4))
 
-n, bins, patches = plt.hist(x, n_bins, normed=1,
-                            histtype='step', cumulative=True)
+# plot the cumulative histograme
+n, bins, patches = ax.hist(x, n_bins, normed=1, histtype='step',
+                           cumulative=True, label='Empirical')
 
 # Add a line showing the expected distribution.
 y = mlab.normpdf(bins, mu, sigma).cumsum()
 y /= y[-1]
-plt.plot(bins, y, 'k--', linewidth=1.5)
+
+ax.plot(bins, y, 'k--', linewidth=1.5, label='Theoretical')
 
 # Overlay a reversed cumulative histogram.
-plt.hist(x, bins=bins, normed=1, histtype='step', cumulative=-1)
+ax.hist(x, bins=bins, normed=1, histtype='step', cumulative=-1,
+        label='Reversed emp.')
 
-plt.grid(True)
-plt.title('cumulative step')
+# tidy up the figure
+ax.grid(True)
+ax.legend(loc='right')
+ax.set_title('Cumulative step histograms')
+ax.set_xlabel('Annual rainfall (mm)')
+ax.set_ylabel('Likelihood of occurrence')
 
 plt.show()
diff --git a/examples/statistics/histogram_demo_features.py b/examples/statistics/histogram_demo_features.py
index 681cb0f7e15f..c662f36ac5e1 100644
--- a/examples/statistics/histogram_demo_features.py
+++ b/examples/statistics/histogram_demo_features.py
@@ -1,35 +1,49 @@
 """
-Demo of the histogram (hist) function with a few features.
+=========================================================
+Demo of the histogram (hist) function with a few features
+=========================================================
 
-In addition to the basic histogram, this demo shows a few optional features:
+In addition to the basic histogram, this demo shows a few optional
+features:
 
     * Setting the number of data bins
-    * The ``normed`` flag, which normalizes bin heights so that the integral of
-      the histogram is 1. The resulting histogram is a probability density.
+    * The ``normed`` flag, which normalizes bin heights so that the
+      integral of the histogram is 1. The resulting histogram is an
+      approximation of the probability density function.
     * Setting the face color of the bars
     * Setting the opacity (alpha value).
 
+Selecting different bin counts and sizes can significantly affect the
+shape of a histogram. The Astropy docs have a great section on how to
+select these parameters:
+http://docs.astropy.org/en/stable/visualization/histogram.html
 """
+
 import numpy as np
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 
+np.random.seed(0)
 
 # example data
 mu = 100  # mean of distribution
 sigma = 15  # standard deviation of distribution
-x = mu + sigma * np.random.randn(10000)
+x = mu + sigma * np.random.randn(437)
 
 num_bins = 50
+
+fig, ax = plt.subplots()
+
 # the histogram of the data
-n, bins, patches = plt.hist(x, num_bins, normed=1)
+n, bins, patches = ax.hist(x, num_bins, normed=1)
+
 # add a 'best fit' line
 y = mlab.normpdf(bins, mu, sigma)
-plt.plot(bins, y, '--')
-plt.xlabel('Smarts')
-plt.ylabel('Probability density')
-plt.title(r'Histogram of IQ: $\mu=100$, $\sigma=15$')
+ax.plot(bins, y, '--')
+ax.set_xlabel('Smarts')
+ax.set_ylabel('Probability density')
+ax.set_title(r'Histogram of IQ: $\mu=100$, $\sigma=15$')
 
 # Tweak spacing to prevent clipping of ylabel
-plt.subplots_adjust(left=0.15)
+fig.tight_layout()
 plt.show()
diff --git a/examples/statistics/histogram_demo_histtypes.py b/examples/statistics/histogram_demo_histtypes.py
index 898b90cae17f..5a6f9fe08a8e 100644
--- a/examples/statistics/histogram_demo_histtypes.py
+++ b/examples/statistics/histogram_demo_histtypes.py
@@ -1,17 +1,25 @@
 """
-Demo of the histogram (hist) function with different ``histtype`` settings.
+================================================================
+Demo of the histogram function's different ``histtype`` settings
+================================================================
 
 * Histogram with step curve that has a color fill.
-* Histogram with with unequal bin widths.
+* Histogram with custom and unequal bin widths.
 
+Selecting different bin counts and sizes can significantly affect the
+shape of a histogram. The Astropy docs have a great section on how to
+select these parameters:
+http://docs.astropy.org/en/stable/visualization/histogram.html
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
+np.random.seed(0)
 
 mu = 200
 sigma = 25
-x = mu + sigma*np.random.randn(10000)
+x = np.random.normal(mu, sigma, size=100)
 
 fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(8, 4))
 
@@ -23,5 +31,5 @@
 ax1.hist(x, bins, normed=1, histtype='bar', rwidth=0.8)
 ax1.set_title('unequal bins')
 
-plt.tight_layout()
+fig.tight_layout()
 plt.show()
diff --git a/examples/statistics/histogram_demo_multihist.py b/examples/statistics/histogram_demo_multihist.py
index 57dd2ac00691..04619cd670e2 100644
--- a/examples/statistics/histogram_demo_multihist.py
+++ b/examples/statistics/histogram_demo_multihist.py
@@ -1,22 +1,31 @@
 """
-Demo of the histogram (hist) function with multiple data sets.
+=============================================================
+Demo of the histogram (hist) function with multiple data sets
+=============================================================
 
 Plot histogram with multiple sample sets and demonstrate:
 
     * Use of legend with multiple sample sets
     * Stacked bars
-    * Step curve with a color fill
+    * Step curve with no fill
     * Data sets of different sample sizes
+
+Selecting different bin counts and sizes can significantly affect the
+shape of a histogram. The Astropy docs have a great section on how to
+select these parameters:
+http://docs.astropy.org/en/stable/visualization/histogram.html
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
+np.random.seed(0)
 
 n_bins = 10
 x = np.random.randn(1000, 3)
 
 fig, axes = plt.subplots(nrows=2, ncols=2)
-ax0, ax1, ax2, ax3 = axes.flat
+ax0, ax1, ax2, ax3 = axes.flatten()
 
 colors = ['red', 'tan', 'lime']
 ax0.hist(x, n_bins, normed=1, histtype='bar', color=colors, label=colors)
@@ -26,13 +35,13 @@
 ax1.hist(x, n_bins, normed=1, histtype='bar', stacked=True)
 ax1.set_title('stacked bar')
 
-ax2.hist(x, n_bins, histtype='step', stacked=True, fill=True)
-ax2.set_title('stepfilled')
+ax2.hist(x, n_bins, histtype='step', stacked=True, fill=False)
+ax2.set_title('stack step (unfilled)')
 
 # Make a multiple-histogram of data-sets with different length.
 x_multi = [np.random.randn(n) for n in [10000, 5000, 2000]]
 ax3.hist(x_multi, n_bins, histtype='bar')
 ax3.set_title('different sample sizes')
 
-plt.tight_layout()
+fig.tight_layout()
 plt.show()
diff --git a/examples/statistics/multiple_histograms_side_by_side.py b/examples/statistics/multiple_histograms_side_by_side.py
index 3b8a1418d06d..b150f8b53197 100644
--- a/examples/statistics/multiple_histograms_side_by_side.py
+++ b/examples/statistics/multiple_histograms_side_by_side.py
@@ -1,15 +1,34 @@
-
 """
+=======================================================
 Demo of how to produce multiple histograms side by side
+=======================================================
+
+This example plots horizontal histograms of different samples along
+a categorical x-axis. Additionally, the histograms are plotted to
+be symmetrical about their x-position, thus making them very similar
+to violin plots.
+
+To make this highly specialized plot, we can't use the standard ``hist``
+method. Instead we use ``barh`` to draw the horizontal bars directly. The
+vertical positions and lengths of the bars are computed via the
+``np.histogram`` function. The histograms for all the samples are
+computed using the same range (min and max values) and number of bins,
+so that the bins for each sample are in the same vertical positions.
+
+Selecting different bin counts and sizes can significantly affect the
+shape of a histogram. The Astropy docs have a great section on how to
+select these parameters:
+http://docs.astropy.org/en/stable/visualization/histogram.html
 """
 
 import numpy as np
 import matplotlib.pyplot as plt
 
+np.random.seed(0)
 number_of_bins = 20
 
 # An example of three data sets to compare
-number_of_data_points = 1000
+number_of_data_points = 387
 labels = ["A", "B", "C"]
 data_sets = [np.random.normal(0, 1, number_of_data_points),
              np.random.normal(6, 1, number_of_data_points),
@@ -17,20 +36,22 @@
 
 # Computed quantities to aid plotting
 hist_range = (np.min(data_sets), np.max(data_sets))
-binned_data_sets = [np.histogram(d, range=hist_range, bins=number_of_bins)[0]
-                    for d in data_sets]
+binned_data_sets = [
+    np.histogram(d, range=hist_range, bins=number_of_bins)[0]
+    for d in data_sets
+]
 binned_maximums = np.max(binned_data_sets, axis=1)
 x_locations = np.arange(0, sum(binned_maximums), np.max(binned_maximums))
 
 # The bin_edges are the same for all of the histograms
 bin_edges = np.linspace(hist_range[0], hist_range[1], number_of_bins + 1)
-centers = .5 * (bin_edges + np.roll(bin_edges, 1))[:-1]
+centers = 0.5 * (bin_edges + np.roll(bin_edges, 1))[:-1]
 heights = np.diff(bin_edges)
 
 # Cycle through and plot each histogram
-ax = plt.subplot(111)
+fig, ax = plt.subplots()
 for x_loc, binned_data in zip(x_locations, binned_data_sets):
-    lefts = x_loc - .5 * binned_data
+    lefts = x_loc - 0.5 * binned_data
     ax.barh(centers, binned_data, height=heights, left=lefts)
 
 ax.set_xticks(x_locations)
diff --git a/examples/statistics/violinplot_demo.py b/examples/statistics/violinplot_demo.py
index 1bfcccf061e0..5aaf8ea7699a 100644
--- a/examples/statistics/violinplot_demo.py
+++ b/examples/statistics/violinplot_demo.py
@@ -1,5 +1,19 @@
 """
-Demo of the new violinplot functionality
+==================================
+Demo of the basics of violin plots
+==================================
+
+Violin plots are similar to histograms and box plots in that they show
+an abstract representation of the probability distribution of the
+sample. Rather than showing counts of data points that fall into bins
+or order statistics, violin plots use kernel density estimation (KDE) to
+compute an empirical distribution of the sample. That computation
+is controlled by several parameters. This example demostrates how to
+modify the number of points at which the KDE is evaluated (``points``)
+and how to modify the band-width of the KDE (``bw_method``).
+
+For more information on violin plots and KDE, the scikit-learn docs
+have a great section: http://scikit-learn.org/stable/modules/density.html
 """
 
 import random
@@ -9,20 +23,20 @@
 # fake data
 fs = 10  # fontsize
 pos = [1, 2, 4, 5, 7, 8]
-data = [np.random.normal(size=100) for i in pos]
+data = [np.random.normal(0, std, size=100) for std in pos]
 
 fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(6, 6))
 
-axes[0, 0].violinplot(data, pos, points=20, widths=0.1,
+axes[0, 0].violinplot(data, pos, points=20, widths=0.3,
                       showmeans=True, showextrema=True, showmedians=True)
 axes[0, 0].set_title('Custom violinplot 1', fontsize=fs)
 
-axes[0, 1].violinplot(data, pos, points=40, widths=0.3,
+axes[0, 1].violinplot(data, pos, points=40, widths=0.5,
                       showmeans=True, showextrema=True, showmedians=True,
                       bw_method='silverman')
 axes[0, 1].set_title('Custom violinplot 2', fontsize=fs)
 
-axes[0, 2].violinplot(data, pos, points=60, widths=0.5, showmeans=True,
+axes[0, 2].violinplot(data, pos, points=60, widths=0.7, showmeans=True,
                       showextrema=True, showmedians=True, bw_method=0.5)
 axes[0, 2].set_title('Custom violinplot 3', fontsize=fs)