def METHOD_NAME(ax, theta):
  dim = theta.shape[1]
  if (dim == 1):
    return
  for i in range(dim):
    for j in range(i):
      # returns bin values, bin edges and bin edges
      H, xe, ye = np.histogram2d(theta[:, j], theta[:, i], 10, density=True)
      # plot and interpolate data
      ax[i, j].imshow(
          H.T,
          aspect="auto",
          interpolation='spline16',
          origin='lower',
          extent=np.hstack((ax[j, j].get_xlim(), ax[i, i].get_xlim())),
          cmap=plt.get_cmap('jet'))
      if i < theta.shape[1] - 1:
        ax[i, j].set_xticklabels([])
      if j > 0:
        ax[i, j].set_yticklabels([])