# F_curve.py

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

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

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

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

p = 0.05
df = 12

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

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

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

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

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

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

ax.set_yticklabels([])
ax.set_ylim([0,0.1])
ax.set_xlim([0,40])
ax.set_title("$\chi^2$ Probability Density Function")

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