code: def average(self, axis):
          
          axis = self.get_axis_number(axis)
          avg_hist = np.ma.average(self.all_axis_bin_centers(axis),
                                   weights=self.histogram, axis=axis)
          if self.dimensions == 2:
              new_hist = Hist1d
          else:
              new_hist = Histdd
          return new_hist.from_histogram(histogram=avg_hist,
                                         bin_edges=itemgetter(*self.other_axes(axis))(self.bin_edges),
                                         axis_names=self.axis_names_without(axis))
docs: Returns d-1 dimensional histogram of (estimated) mean value of axis
          NB this is very different from averaging over the axis!!!