Why simulated and calculated t-distribution pdf at degree-of-freedom=1 doesn't match See the attached graph, I am simulating the pdf of t-distribution at different degrees of freedom (dof), when dof is low (e.g. dof=1), why don't they match? Is it because my experiment wasn't done correctly, or is that expected?
The code
def t_dist_pdf(t, dof):
    p1 = special.gamma((dof + 1) / 2)
    p2 = np.sqrt(dof * np.pi) * special.gamma(dof / 2)
    p3 = (1 + t ** 2 / dof) ** (- (dof + 1) / 2)
    return p1 / p2 * p3

fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(12, 6))
axes = axes.ravel()

bins = np.linspace(-5, 5, 100)
num = int(1e6)

for k, dof in enumerate(range(1, 7)):
    ax = axes[k]
    t_dist_data = np.random.standard_t(df=dof, size=num)
    ax.hist(t_dist_data, bins=bins, normed=True, histtype='step', label='sim. $t$-dist pdf')

    ax.plot(xs, t_dist_pdf(xs, dof), lw=1, alpha=1, label='$t$-dist pdf')
    ax.plot(xs, stats.t.pdf(xs, dof), lw=5, alpha=0.2, label='scipy $t$-dist pdf')
    ax.plot(xs, stats.norm.pdf(xs), lw=1, alpha=1, label='normal pdf')

    ax.legend(fontsize=8)
    ax.set_title('$t$-distribution dof={0}'.format(dof))
plt.tight_layout()

I verified my implementation of t-dist pdf by comparing it to scipy's.
The result. 
The X and Y axes of all subplots are the same.
PDF:

Update: CDF:

Version info:
python: 3.5.3,
numpy: 1.12.1,
scipy: 0.19.1
 A: wuber pointed out the problem,

In particular, by limiting it to the range [−5,5][−5,5] and requiring
  the plot to integrate to unity, you are estimating a truncated
  t-distribution. For the smaller values of df, an appreciable
  proportion of its probability lies beyond that range (12.6% for df=1,
  3.8% for df=2, etc). Accordingly, your plot of the simulated distribution for df=1 should be about 1/(1-12.6%) = 14% too high, etc.
  Equivalently, you should renormalize the theoretical expressions to
  account for the truncation.

Here is the fixed code with proper normalization for the theoretical t-distribution PDF before comparing it to the simulated one:
fig, axes = plt.subplots(4, 3, sharex=False, sharey=False, figsize=(12, 15))
axes = axes.ravel()

bins = np.linspace(-5, 5, 100)
delta_bin = bins[1] - bins[0]

xs = (bins[1:] + bins[:-1]) / 2
num = int(1e6)

for k, dof in enumerate(range(1, 7)):
    ax = axes[k]

    # IMPORTANT: used to correct the theorectical pdf 
    # for truncated regions of t dist and gaussian/normal dist
    t_renorm = 1 - stats.t.cdf(bins[0], df=dof) * 2
    g_renorm = 1 - stats.norm.cdf(bins[0]) * 2

    t_dist_data = np.random.standard_t(df=dof, size=num)

    ax.plot(xs, stats.t.pdf(xs, dof) / t_renorm, lw=5, alpha=0.2, color='cyan', label='scipy pdf')
    ax.plot(xs, t_dist_pdf(xs, dof) / t_renorm, lw=1, alpha=1, color='red', label='self implemented pdf')
    ax.plot(xs, stats.norm.pdf(xs) / g_renorm, lw=1, alpha=1, color='black', label='normal pdf')
    ax.hist(t_dist_data, bins=bins, normed=True, histtype='step', color='blue', label='simulated pdf')

    ax.legend(fontsize=10)
    ax.set_title('$t$-distribution PDF dof={0}'.format(dof))
    ax.set_xlim(bins[0], bins[-1])
    ax.set_ylim(0, 0.6)


    ax2 = axes[k + 6]
    ax2.plot(xs, cdf(stats.t.pdf(xs, dof) / t_renorm, delta_bin), lw=5, alpha=0.2, color='cyan', label='scipy cdf')
    ax2.plot(xs, cdf(t_dist_pdf(xs, dof) / t_renorm, delta_bin), lw=1, alpha=1, color='black', label='simulated cdf')
    ax2.plot(xs, cdf(stats.norm.pdf(xs) / g_renorm, delta_bin), lw=1, alpha=1, color='black', label='normal cdf')
    ax2.legend(loc='lower right', fontsize=10)
    ax2.set_title('$t$-distribution CDF dof={0}'.format(dof))
    ax2.set_xlim(bins[0], bins[-1])
    ax2.set_ylim(-0.05, 1.05)
plt.tight_layout()

The output

Now, everything matches.
