How to use the corner.hist2d function in corner

def _run_hist2d(nm, N=50000, seed=1234, **kwargs):
    print(" .. {0}".format(nm))

    if not os.path.exists(FIGURE_PATH):
        os.makedirs(FIGURE_PATH)

    # Generate some fake data.
    np.random.seed(seed)
    x = np.random.randn(N)
    y = np.random.randn(N)

    fig, ax = pl.subplots(1, 1, figsize=(8, 8))
    corner.hist2d(x, y, ax=ax, **kwargs)
    fig.savefig(os.path.join(FIGURE_PATH, "hist2d_{0}.png".format(nm)))
    pl.close(fig)

def _run_hist2d(nm, N=50000, seed=1234, **kwargs):
    # Generate some fake data.
    np.random.seed(seed)
    x = np.random.randn(N)
    y = np.random.randn(N)

    fig, ax = pl.subplots(1, 1, figsize=(8, 8))
    corner.hist2d(x, y, ax=ax, **kwargs)

'''.split()))
plt.rc('font', family='serif')

f = pyfits.open("../make_lamost_catalog/table_for_paper.fits")
a = f[1].data
f.close()

feh = a['cannon_m_h']
am = a['cannon_alpha_m']
snr = a['snrg']

choose = snr > 50

fig, axarr = plt.subplots(2,1, sharex=True, sharey=True, figsize=(8,10))
hist2d(feh[choose], am[choose], ax=axarr[0], bins=100, range=[[-2.2,.9],[-0.2,0.5]])
hist2d(feh, am, ax=axarr[1], bins=100, range=[[-2.2,.9],[-0.2,0.5]])
axarr[1].set_xlabel("[Fe/H] (dex)" + " from Cannon/LAMOST", fontsize=16)
#axarr[1].set_ylabel(r"$\mathrm{[\alphaup/M]}$" + " (dex) from Cannon/LAMOST", fontsize=16)
fig.text(0.04, 0.5, r"$\mathrm{[\alphaup/M]}$" + " (dex) from Cannon/LAMOST", fontsize=16,
        va = 'center', rotation='vertical')
labels = [r"Objects with SNR \textgreater 50", r"All Objects"]
props = dict(boxstyle='round', facecolor='white')
for i,ax in enumerate(axarr):
    ax.tick_params(axis='x', labelsize=16)
    ax.tick_params(axis='y', labelsize=16)
    ax.text(0.05, 0.85, labels[i], 
            horizontalalignment='left', verticalalignment='bottom', transform=ax.transAxes,
            fontsize=16, bbox=props)
#plt.show()
plt.savefig("feh_alpha.png")

else:
                axis.set_ylabel(params[i])

            if i == j:                    
                axis.hist(fit1.posterior["samples"][:,param_indices1[i]], bins=20, color=colour1_2, normed=True, histtype="stepfilled", edgecolor=colour1, range=[np.min([lims1[i,0], lims2[i,0]]), np.max([lims1[i,1], lims2[i,1]])], zorder=6, alpha=alpha1, lw=2)
                axis.hist(fit2.posterior["samples"][:,param_indices2[i]], bins=20, color=colour2_2, normed=True, histtype="stepfilled", edgecolor=colour2, range=[np.min([lims1[i,0], lims2[i,0]]), np.max([lims1[i,1], lims2[i,1]])], zorder=7, alpha=alpha2, lw=2)

                axis.set_xlim(np.min([lims1[i,0], lims2[i,0]]), np.max([lims1[i,1], lims2[i,1]]))

                if truths is not None:
                    axis.axvline(truths[i], lw=2, color="#4682b4", zorder=8)

            else:
                corner.hist2d(fit1.posterior["samples"][:,param_indices1[j]], fit1.posterior["samples"][:,param_indices1[i]], ax=axis, smooth=1.5, range=(lims1[j,:], lims1[i,:]), color=colour1)#, plot_density=False)
                corner.hist2d(fit2.posterior["samples"][:,param_indices2[j]], fit2.posterior["samples"][:,param_indices2[i]], ax=axis, smooth=1.5, range=(lims2[j,:], lims2[i,:]), color=colour2)#, plot_density=False)

                axis.set_xlim(np.min([lims1[j,0], lims2[j,0]]), np.max([lims1[j,1], lims2[j,1]]))
                axis.set_ylim(np.min([lims1[i,0], lims2[i,0]]), np.max([lims1[i,1], lims2[i,1]]))

                if truths is not None:
                    axis.axhline(truths[i], lw=2, color="#4682b4", zorder=8)
                    axis.axvline(truths[j], lw=2, color="#4682b4", zorder=8)
                    axis.plot(truths[j], truths[i], "s", color="#4682b4", zorder=8)


        
    sfh_ax = fig.add_axes([0.65, 0.59, 0.32, 0.15], zorder=10)
    sfr_ax = fig.add_axes([0.82, 0.82, 0.15, 0.15], zorder=10)
    tmw_ax = fig.add_axes([0.65, 0.82, 0.15, 0.15], zorder=10)

system.compute_orbits(time)
            duration_zoom[n,i] = \
                dt * len(np.where(system.bodies[1].occultor == 4)[0])
    
    # Plot the heatmap
    fig = pl.figure(figsize = (6, 8))
    ax = [pl.subplot2grid((3,1), (0,0), rowspan = 2),
                pl.subplot2grid((3,1), (2,0))]

    # Top plot
    duration = duration.reshape(-1)
    Omega = np.concatenate([Omega for i in range(niter)])
    iszero = np.where(duration == 0)
    duration = np.delete(duration, iszero)
    Omega = np.delete(Omega, iszero)
    corner.hist2d(duration, Omega, bins = 100, ax = ax[0])

    # Zoom inset
    duration_zoom = duration_zoom.reshape(-1)
    Omega_zoom = np.concatenate([Omega_zoom for i in range(niter)])
    iszero = np.where(duration_zoom == 0)
    duration_zoom = np.delete(duration_zoom, iszero)
    Omega_zoom = np.delete(Omega_zoom, iszero)
    corner.hist2d(duration_zoom, Omega_zoom, bins = 100, ax = ax[1])

    # Appearance
    ax[0].set_ylabel(r'$\mathbf{\Delta\Omega}$ [degrees]', 
                     fontsize = 14, fontweight = 'bold')
    ax[1].set_ylabel(r'$\mathbf{\Delta\Omega}$ [degrees]', 
                     fontsize = 14, fontweight = 'bold')
    ax[1].set_xlabel('Duration of mutual transit [minutes]', 
                     fontsize = 14, fontweight = 'bold')

n, _, _ = axt[p].hist(y, bins = bins, orientation="vertical", histtype='step', 
                              fill=False, color = 'k', lw = 1, weights = weights)
        ax[p].set_ylim(0, 1.2 * maxn2) 
        axt[p].set_ylim(maxn1 - 1.2 * maxn2 / 8, maxn1 + 1.2 * maxn2 / 8)
        axt[p].set_yticks([maxn1])
        axt[p].set_yticklabels(['%.2f' % maxn1])
        ax[p].yaxis.set_label_coords(-0.15, 0.67)
      else:
        ax[p].set_ylim(0, 1.2 * maxn1)
      
    # Corner
    blobs = blobs[nburn:,:,np.array([-2,-1])].reshape(nwalk * (nsteps - nburn), 2)
    water = blobs[:,0]
    oxygen = blobs[:,1]
    matplotlib.rcParams['lines.linewidth'] = 1
    corner.hist2d(water, oxygen, ax = ax[2], plot_density = False, plot_contours = False, data_kwargs = {'alpha': 0.002})
    
    # Appearance
    if label is not None:
      ax[0].annotate(label, xy = (-0.3, 0.5), xycoords = 'axes fraction', ha = 'right', va = 'center', fontsize = 22)
    ax[0].set_xlim(wxmin, wxmax)
    ax[0].get_xaxis().set_major_locator(MaxNLocator(6))
    ax[1].set_xlim(oxmin, oxmax)
    ax[1].get_xaxis().set_major_locator(MaxNLocator(5))
    ax[2].set_xticks(ax[0].get_xticks())
    ax[2].set_xlim(wxmin, wxmax)
    ax[2].set_yticks(ax[1].get_xticks())
    ax[2].set_ylim(oxmin, oxmax)
    ax[2].set_xlabel('Water (TO)', fontsize = 14)
    ax[2].set_ylabel('Oxygen (bar)', fontsize = 14)
    for n in [0,1]:
      ax[n].yaxis.set_major_formatter(FuncFormatter(lambda x, p : '%.2f' % x))

# Top plot
    duration = duration.reshape(-1)
    Omega = np.concatenate([Omega for i in range(niter)])
    iszero = np.where(duration == 0)
    duration = np.delete(duration, iszero)
    Omega = np.delete(Omega, iszero)
    corner.hist2d(duration, Omega, bins = 100, ax = ax[0])

    # Zoom inset
    duration_zoom = duration_zoom.reshape(-1)
    Omega_zoom = np.concatenate([Omega_zoom for i in range(niter)])
    iszero = np.where(duration_zoom == 0)
    duration_zoom = np.delete(duration_zoom, iszero)
    Omega_zoom = np.delete(Omega_zoom, iszero)
    corner.hist2d(duration_zoom, Omega_zoom, bins = 100, ax = ax[1])

    # Appearance
    ax[0].set_ylabel(r'$\mathbf{\Delta\Omega}$ [degrees]', 
                     fontsize = 14, fontweight = 'bold')
    ax[1].set_ylabel(r'$\mathbf{\Delta\Omega}$ [degrees]', 
                     fontsize = 14, fontweight = 'bold')
    ax[1].set_xlabel('Duration of mutual transit [minutes]', 
                     fontsize = 14, fontweight = 'bold')
    ax[0].set_ylim(0, 180)
    ax[1].set_ylim(0, 3)

    return fig, ax

plt.rc('font', family='serif')

direc = '/users/annaho/Data/LAMOST/Mass_And_Age'
f = pyfits.open("%s/Ho2016b_Catalog.fits" %direc)
a = f[1].data
f.close()

feh = a['MH']
am = a['AM']
snr = a['SNR']
chisq = a['RedChisq']

choose = snr > 20

fig, ax = plt.subplots(1,1, sharex=True, sharey=True, figsize=(8,5))
hist2d(feh[choose], am[choose], ax=ax, bins=100, range=[[-2.2,.9],[-0.2,0.5]])
ax.set_xlabel("[Fe/H] (dex)" + " from Cannon/LAMOST", fontsize=16)
fig.text(
        0.04, 0.5, r"$\mathrm{[\alpha/M]}$" + " (dex) from Cannon/LAMOST", 
        fontsize=16, va = 'center', rotation='vertical')
label = r"Objects with SNR \textgreater 20"
props = dict(boxstyle='round', facecolor='white')
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax.text(0.05, 0.85, label, 
        horizontalalignment='left', verticalalignment='bottom', 
        transform=ax.transAxes, fontsize=16, bbox=props)
#plt.show()
plt.savefig("feh_alpha.png")