import numpy as np
from scipy import stats as st, optimize as opt
from matplotlib import pyplot as plt
%matplotlib inline

normal_prob = lambda x: st.norm.cdf(x) - st.norm.cdf(-x)
normal_range = lambda x: (st.norm.ppf((1-x)/2), st.norm.isf((1-x)/2))

print('The probability for two sigmas is:', normal_prob(2))
print('The central range for 95% probability is:', normal_range(0.95))

The probability for two sigmas is: 0.954499736104
The central range for 95% probability is: (-1.959963984540054, 1.959963984540054)

def chi2_prob(x, df):
    std = np.sqrt(2*df)
    return st.chi2.cdf(df+std*x, df) - st.chi2.cdf(df-std*x, df)

print('The probability for two sigmas is:', chi2_prob(2, 50))
print('The difference with the normal aproximation is:', chi2_prob(2, 50)-normal_prob(2))

The probability for two sigmas is: 0.956461109955
The difference with the normal aproximation is: 0.00196137385127

res = opt.minimize_scalar(lambda x: np.abs(chi2_prob(x, 50)-0.95))

print('El intervalo de 95% de probabilidad con el mismo ancho para ambos lados es:', res.x,
      'sigmas (para la normal era', normal_range(0.95)[1], 'sigmas).')

print('Es decir, el intervalo es:', (50-10*res.x, 50+10*res.x))
print('mientras que en la aproximación normal es:', (50-10*normal_range(0.95)[1], 50+10*normal_range(0.95)[1]))

El intervalo de 95% de probabilidad con el mismo ancho para ambos lados es: 1.93922366658 sigmas (para la normal era 1.95996398454 sigmas).
Es decir, el intervalo es: (30.607763334179761, 69.392236665820235)
mientras que en la aproximación normal es: (30.400360154599461, 69.599639845400532)

sigma = np.linspace(0, 6, 100)
plt.plot(sigma, chi2_prob(sigma, 50) - normal_prob(sigma))
plt.vlines(5, -0.002, 0.006, linestyle='--')
plt.ylim(-0.002, 0.006)
plt.xlabel('Std. dev. from mean')
plt.ylabel('Difference of probabilities')

print('Nota: A partir de los 5 sigmas, la chi^2 de 50 grados de libertad es nula (del lado izquierdo) ya que estamos en números negativos.')

Nota: A partir de los 5 sigmas, la chi^2 de 50 grados de libertad es nula (del lado izquierdo) ya que estamos en números negativos.

for sigma in np.linspace(0, 5, 20):
    plt.plot([central_probability_chi(sigma, df) - central_probability(sigma) for df in range(1, 100)])

plt.ylim(-0.05, 0.05)

(-0.05, 0.05)

x = np.linspace(0, 150, 1000)

f = st.chi2(df=50)

print('Probabilidad de que salga 80 o más:', f.sf(80))

plt.plot(x, f.pdf(x))
plt.fill_between(x[x>80], f.pdf(x[x>80]), color='r')

Probabilidad de que salga 80 o más: 0.00448265656557

<matplotlib.collections.PolyCollection at 0x7fcc27d4c828>