code: def get_random(self, size=10):
          
          # Sample random bin centers
          bin_centers_ravel = np.array(np.meshgrid(*self.bin_centers(),
                                                   indexing='ij')).reshape(self.dimensions, -1).T
          hist_ravel = self.histogram.ravel()
          hist_ravel = hist_ravel.astype(np.float) / np.nansum(hist_ravel)
          result = bin_centers_ravel[np.random.choice(len(bin_centers_ravel),
                                                      p=hist_ravel,
                                                      size=size)]
  
          # Randomize the position inside the bin
          for dim_i in range(self.dimensions):
              bin_edges = self.bin_edges[dim_i]
              bin_widths = np.diff(bin_edges)
  
              # Note the - 1: for the first bin's bin center, searchsorted gives 1, but we want 0 here:
              index_of_bin = np.searchsorted(bin_edges, result[:, dim_i]) - 1
              result[:, dim_i] += (np.random.rand(size) - 0.5) * bin_widths[index_of_bin]
  
          return result
docs: Returns (size, n_dim) array of random variates from the histogram.
          Inside the bins, a uniform distribution is assumed
          Note this assumes the histogram is an 'events per bin', not a pdf.
          TODO: test more.