diff --git a/fluents/plots.py b/fluents/plots.py
index 617a1f8..620306b 100644
--- a/fluents/plots.py
+++ b/fluents/plots.py
@@ -575,7 +575,8 @@ class HistogramPlot(Plot):
 
         # Set default paramteters
         if not kw.has_key('bins'):
-            kw['bins'] = 20
+            kw['bins'] = self._get_binsize()
+            print kw['bins']
 
         # Initial draw
         self.axes.grid(False)
@@ -598,7 +599,7 @@ class HistogramPlot(Plot):
             self._toolbar.freeze_button.set_sensitive(False)
             self._toolbar.set_mode_sensitive('select', False)
             self._toolbar.set_mode_sensitive('lassoselect', False)
-    
+
     def rectangle_select_callback(self, x1, y1, x2, y2, key):
         if self.current_dim == None: return
         # make (x1, y1) the lower left corner
@@ -653,6 +654,16 @@ class HistogramPlot(Plot):
                 patch.set_facecolor((r,g,b,1))
         self.canvas.draw()
 
+    def _get_binsize(self, min_bins=2, max_bins=100):
+        """ Automatic bin selection, as described by Shimazaki."""
+        bin_vec = scipy.arange(min_bins, max_bins, 1)
+        D = self._data.ptp()/bin_vec
+        cost = scipy.empty((bin_vec.shape[0],), 'f')
+        for i, bins in enumerate(bin_vec):
+            count, lims = scipy.histogram(self._data, bins)
+            cost[i] = (2*count.mean() - count.var())/(D[i]**2)
+        best_bin_size = bin_vec[scipy.argmin(cost)]
+        return best_bin_size
 
 class BarPlot(Plot):
     """Bar plot.