# t_curve.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t

#########################################################################################

def indices(a, func):
   return [i for (i, val) in enumerate(a) if func(val)]

#########################################################################################

p = 0.05
df = 12

value = t.ppf(1-p/2, df)

x = np.arange(-value, value, 0.01)

x_all = np.arange(-10, 10, 0.001) # entire range of x, both in and out of spec

y = t.pdf(x, df=12)
y2 = t.pdf(x_all,df=12)

# build the plot
fig, ax = plt.subplots(figsize=(9,6))
# plt.style.use('fivethirtyeight')
ax.plot(x_all, y2, color='b')

j = indices(x_all,  lambda y: y < -value)
k = indices(x_all,  lambda y: y > value)
ax.fill_between(x_all[j], t.pdf(x_all[j], df), 0, alpha=0.3, color='b')
ax.fill_between(x_all[k], t.pdf(x_all[k], df), 0, alpha=0.3, color='b')
ax.fill_between(x_all, y2, 0, color='b', alpha=0.1)
ax.set_xlim([-4,4])
ax.set_xlabel('Number of Standard Deviations from Mean 0')
ax.set_ylim([0, 0.4])
ax.set_yticklabels([])
ax.set_title("Student's $t$ Probability Density Function")

plt.savefig('t_curve.png', dpi=200, bbox_inches='tight', transparent=True)
plt.show()