# F_curve.py

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

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

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

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

p = 0.05
df = 29

value1 = f.ppf(p/2, df, df)
value2 = f.ppf(1-p/2, df, df)

x = np.arange(value1, value2, 0.01)
x_all = np.arange(0.0001, 3, 0.001) # entire range of x, both in and out of spec

y = f.pdf(x, df, df)
y2 = f.pdf(x_all, df, df)

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

j = indices(x_all,  lambda y: y < value1)
k = indices(x_all,  lambda y: y > value2)

ax.fill_between(x_all[j], f.pdf(x_all[j], df, df), 0, alpha=0.3, color='g')
ax.fill_between(x_all[k], f.pdf(x_all[k], df, df), 0, alpha=0.3, color='g')
ax.fill_between(x_all, y2, 0, alpha=0.1, color='g')

ax.set_xlim([0,3])
ax.set_ylim([0,1.2])
ax.set_yticklabels([])
ax.set_title("Snedecor's  $F$ Probability Density Function")

plt.savefig('f_curve.png', dpi=200, bbox_inches='tight', transparent=True, \
             facecolor='w', edgecolor='none')
plt.show()