#!/usr/bin/env python
# Get python from http://www.python.org/

# Example code to implement Gaussian and Student's t-distributions

# Standard python library modules
import math
import sys

# Get scipy from http://www.scipy.org/
from scipy import special

class Gauss(object):
    def pdf(self, z):
        return 1.0/math.sqrt(2.0*math.pi) * math.exp(-z**2/2)
    
    def cdf(self, z):
        return 0.5*(special.erfc(-z/math.sqrt(2.0)))
        
    def cdf_inv(self, p):
        if p <= 0 or p >= 1:
            raise ValueError("p must be between 0 and 1")
        c0 = 2.515517
        c1 = 0.802853
        c2 = 0.010328
        d1 = 1.432788
        d2 = 0.189269
        d3 = 0.001308

        sign = -1.
        if (p > 0.5):
                sign = 1.
                p = 1. - p
        arg = -2.* math.log(p)
        t = math.sqrt(arg)
        g = t - (c0 + t*(c1 + t*c2)) / (1. + t*(d1 + t*(d2 + t*d3)))
        return sign*g

class T(object):
    def __init__(self, nu):
        self._nu = float(nu)
        self.calcPrefactor()

    def calcPrefactor(self):
        self._n12 = (self._nu+1.0)/2.0
        self._g12 = special.gamma(self._n12)
        self._gn2 = special.gamma(self._nu/2.0)
        t = self._g12/self._gn2
        self._spn = math.sqrt(self._nu*math.pi)
        self._pdf0 = t/self._spn
        
    def pdf(self, x):
        t = (1.0 + x**2/self._nu) ** (-self._n12)
        return self._pdf0 * t

    def cdf(self, x):
        t = special.hyp2f1(0.5, self._n12, 1.5, -x**2/self._nu)
        u = self._spn * self._gn2
        return 0.5 + x*self._g12 * t/u

    def cdf_inv(self, y):
        if y <= 0 or y >= 1:
            raise ValueError("y must be between 0 and 1")
        if y<0.5:
            xhigh = 0.0
            xlow = -1.0
            while self.cdf(xlow) > y:
                xhigh = xlow
                xlow *= 2
        else:
            xlow = 0.0
            xhigh = 1.0
            while self.cdf(xhigh) < y:
                xlow = xhigh
                xhigh *= 2
        while xhigh - xlow > 1e-6:
            xmid = (xhigh + xlow)/2.0
            if self.cdf(xmid) < y:
                xlow = xmid
            else:
                xhigh = xmid
        return (xhigh + xlow)/2.0

def graph():
    g=Gauss()
    t10=T(10)
    t3=T(3)
    for i in range(0,1000):
        x=(i-500)/100.
        print "%.3f %.7f %.7f %.7f"% (x,g.pdf(x),t10.pdf(x),t3.pdf(x))