def METHOD_NAME(ax, theta):
  dim = theta.shape[1]
  num_bins = 50
  for i in range(dim):
    if (dim == 1):
      ax_loc = ax
    else:
      ax_loc = ax[i, i]
    hist, bins, _ = ax_loc.hist(
        theta[:, i], num_bins, density=True, color='lightgreen', ec='black')
    if i == 0:
      # Rescale hist to scale of theta -> get correct axis titles
      widths = np.diff(bins)
      if (dim > 1):
        hist = hist / np.max(hist) * (
            ax_loc.get_xlim()[1] - ax_loc.get_xlim()[0])
        bottom = ax_loc.get_xlim()[0]
        ax_loc.cla()
        ax_loc.bar(
            bins[:-1],
            hist,
            widths,
            color='lightgreen',
            ec='black',
            bottom=bottom)
        ax_loc.set_ylim(ax_loc.get_xlim())
        ax_loc.set_xticklabels([])
      else:
        ax_loc.cla()
        ax_loc.bar(bins[:-1], hist, widths, color='lightgreen', ec='black')
    elif i == theta.shape[1] - 1:
      ax_loc.set_yticklabels([])
    else:
      ax_loc.set_xticklabels([])
      ax_loc.set_yticklabels([])
    ax_loc.tick_params(axis='both', which='both', length=0)