import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sp
import scipy.integrate as integral
 
col = ['orange', 'purple', 'deepskyblue']
 
X = np.arange(-5, 5, 0.005)
 
plt.clf()
plt.figure(figsize=(4,3.2))
plt.axes([0.17,0.13,0.79,0.8])
plt.hold(True)
 
A = []

def PDF(x,nu):
	Y = -(nu+1)*np.log(1+x**2/nu)/2
	Y += sp.gammaln((nu+1)/2.0)
	Y -= sp.gammaln(nu/2.0)
	Y -= 0.5*np.log(nu*np.pi)
	Y = np.exp(Y)
	return Y

for k,nu in enumerate([1,2,5]):
	Y=[]
	for x in X:
		Y.append(integral.quad(lambda t: PDF(t,nu),-1000,x))
	a = plt.plot(X, Y, '-', color=col[k], lw=1.5)
	A.append(a)

# Degenerate into Normal distribution
from scipy.stats import norm
Y = norm.cdf(X)
a = plt.plot(X, Y, '-', color='black', lw=1.5)
A.append(a)
 
plt.xlabel("x")
plt.ylabel(u"P(X\N{Less-THAN OR EQUAL TO}x)")

bx = plt.legend(A, (r"$\nu=1$", r"$\nu=2$", r"$\nu=5$", r"$\nu=+\infty$"),\
                numpoints=1, handlelen=0.05, handletextpad=0.4,\
                loc="lower right")
bx.draw_frame(False)
plt.xlim(-5,5)
 
plt.savefig("student_t_cdf.pdf")
plt.savefig("student_t_cdf.eps")
plt.savefig("student_t_cdf.svg")