import matplotlib.pyplot as plt
import numpy as np
import matplotlib.gridspec as gridspec
import scipy.special
from scipy.optimize import curve_fit
import seaborn as sns
set(style="ticks", font_scale=1.5) sns.
Fun with histograms
1d and 2d histograms
Introduction
This code produces the figure above. I tried to showcase a few things one can do with 1d and 2d histograms.
The code
= plt.figure(1, figsize=(8,6)) # figsize accepts only inches.
fig
##################################
# Panels on the left of the figure
##################################
= gridspec.GridSpec(2, 2, width_ratios=[1, 0.2], height_ratios=[0.2, 1])
gs =0.10, right=0.50, top=0.95, bottom=0.13,
gs.update(left=0.02, wspace=0.02)
hspace
= 1.0 # standard deviation (spread)
sigma = 0.0 # mean (center) of the distribution
mu = np.random.normal(loc=mu, scale=sigma, size=5000)
x = 2.0 # shape
k = 1.0 # scale
theta = np.random.gamma(shape=k, scale=theta, size=5000)
y
# bottom left panel
= plt.subplot(gs[1, 0])
ax10 = ax10.hist2d(x, y, bins=40, cmap="YlOrRd",
counts, xedges, yedges, image =True)
density= xedges[1] - xedges[0]
dx = yedges[1] - yedges[0]
dy = xedges[:-1] + dx / 2
xvec = yedges[:-1] + dy / 2
yvec r"$x$")
ax10.set_xlabel(r"$y$", rotation="horizontal")
ax10.set_ylabel(-2, 8, r"$p(x,y)$")
ax10.text(min(), xedges.max()])
ax10.set_xlim([xedges.min(), yedges.max()])
ax10.set_ylim([yedges.
# top left panel
= plt.subplot(gs[0, 0])
ax00 = (1.0 / np.sqrt(2.0 * np.pi * sigma ** 2)) * \
gaussian -((xvec - mu) ** 2) / (2.0 * sigma ** 2))
np.exp(= counts.sum(axis=1) * dy
xdist =dx, fill=False,
ax00.bar(xvec, xdist, width='black', alpha=0.8)
edgecolor='black')
ax00.plot(xvec, gaussian, colormin(), xedges.max()])
ax00.set_xlim([xedges.
ax00.set_xticklabels([])
ax00.set_yticks([])"Normal distribution", fontsize=16)
ax00.set_xlabel("top")
ax00.xaxis.set_label_position(r"$p(x)$", rotation="horizontal", labelpad=20)
ax00.set_ylabel(
# bottom right panel
= plt.subplot(gs[1, 1])
ax11 = yvec ** (k - 1.0) * np.exp(-yvec / theta) / \
gamma_dist ** k * scipy.special.gamma(k))
(theta = counts.sum(axis=0) * dx
ydist =dy, fill=False,
ax11.barh(yvec, ydist, height='black', alpha=0.8)
edgecolor='black')
ax11.plot(gamma_dist, yvec, colormin(), yedges.max()])
ax11.set_ylim([yedges.
ax11.set_xticks([])
ax11.set_yticklabels([])"Gamma distribution", fontsize=16)
ax11.set_ylabel("right")
ax11.yaxis.set_label_position(r"$p(y)$")
ax11.set_xlabel("top")
ax11.xaxis.set_label_position(
###################################
# Panels on the right of the figure
###################################
= gridspec.GridSpec(2, 1, width_ratios=[1], height_ratios=[1, 1])
gs2 =0.60, right=0.98, top=0.95, bottom=0.13,
gs2.update(left=0.02, wspace=0.05)
hspace
= np.random.normal(loc=0, scale=1, size=1000)
x = np.random.gamma(shape=2, size=1000)
y
= plt.subplot(gs2[1, 0])
bx10 = plt.subplot(gs2[0, 0])
bx00
= 100
N = np.random.gamma(shape=5, size=N)
a = np.arange(0,15,1.5)
my_bins = bx00.hist(a, bins=my_bins, density=True,
n1, bins1, patches1 ='stepfilled', alpha=0.2, hatch='/')
histtype0, 15])
bx00.set_xlim([0, 0.28])
bx00.set_ylim([
bx00.set_xticklabels([])"plt.hist", fontfamily="monospace")
bx00.set_xlabel("top")
bx00.xaxis.set_label_position(
# the following way is equivalent to plt.hist, but it gives
# the user more flexibility when plotting and analysing the results
= np.histogram(a, bins=my_bins, density=True)
n2, bins2 = bins2[1] - bins2[0]
wid = bx10.plot(bins2[:-1]+wid/2, n2, marker='o', color='red')
red, -1], n2, width=wid, fill=False,
bx10.bar(bins2[:='black', linewidth=3, alpha=0.8, align="edge")
edgecolor0, 15])
bx10.set_xlim([0, 0.28])
bx10.set_ylim([r"np.histogram; plt.bar", fontfamily="monospace")
bx10.set_xlabel(
# best fit
= my_bins[:-1] + wid/2
xdata = n2
ydata def func(x, p1, p2):
return x ** (p1 - 1.0) * np.exp(-x / p2) / (p2 ** p1 * scipy.special.gamma(p1))
= curve_fit(func, xdata, ydata, p0=(1.5, 1.5)) # p0 = initial guess
popt, pcov = popt
p1, p2 = ((ydata - ydata.mean()) ** 2).sum()
SStot = ((ydata - func(xdata, p1, p2)) ** 1).sum()
SSres = 1 - SSres / SStot
Rsquared = np.linspace(0,15,101)
h ='blue', linewidth=2)
bx00.plot(h, func(h, p1, p2), color# dummy plot, just so we can have a legend on the bottom panel
= ax10.plot([100],[100], color='blue', linewidth=2, label="Best fit")
blue, r'Data',r'Best fit, $r^2=${:.2f}'.format(Rsquared)],
bx10.legend([red,blue],[='upper right', frameon=False, handlelength=4,
loc=False, numpoints=3, fontsize=14)
markerfirst
"histograms.png",dpi=300) fig.savefig(