import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.special import erf
from scipy.stats import chi2

Error function and cumulative of the gaussian

$erf(x) = \frac{1}{2\pi} \int_0^x e^{-t^2}dt$

Cumulative of "unit" gaussian $G(\mu=0, \sigma=1)$ is $0.5(1+errf(x))$

# Error function is erf(x) = (1/2pi)*Integral_from_0_to_x of e^{-t^2}dt
# The cumulative of the 'unit' gaussian is = 0.5 * (1+erf(x)) 
x  = np. linspace(-3,3,2000)
ax = plt.subplot(111)
ax.plot(x, 0.5*(1+erf(x)))
ax.grid()
ax.set_xlabel("$\sigma$ from mean")
_ = ax.set_ylim(bottom=0)

---

$\chi^2 = \sum_i^N \frac{(x_i-\mu_i)^2}{\sigma_i^2}$

Number of degrees of freedom = $n = N$
or $n = N-p$ if $\mu_i$ are fit to $p$ parameters
If $x_i$ are $G(\mu_i, \sigma_i)$ then $\chi^2$ follows this distribution

$f(x) dx = \frac{1}{\Gamma(n/2)} x^{n/2-1} e^{-n/2}$

The expectation value (mean) is $<\chi^2> = n$
$\chi^2/n \approx 1$ is a rule-of-thumb for a good fit
But it there is much more to it in the tails.

# Plot chisquare per ndof for different ndof
ndof = [1,2,4,8,15,20]
x  = np. linspace(0,5,1000)
figure = plt.figure(figsize=(10, 4.5))
ax = plt.subplot(1,2,1)
for n in ndof:
    ax.plot(x, chi2.pdf(n*x, df=n), label="ndof ="+str(n))
ax.grid()
ax.set_ylim(0.,0.5)
ax.legend()
ax.set_xlabel("$\chi^2$/ndof")
ax.set_xlim(x[0],x[-1])

ax = plt.subplot(1,2,2)
for n in ndof:
    ax.plot(x, chi2.pdf(n*x, df=n), label="ndof ="+str(n))
ax.grid()
ax.legend()
ax.set_xlabel("$\chi^2$/ndof")
ax.set_yscale('log')
ax.set_ylim(1e-6,1)
_ = ax.set_xlim(x[0],x[-1])