from scipy import *
from scipy import random 
from scipy import stats
import pylab

# Task 1

def sampleExpo(lam, n):
    u = random.random_sample(n)
    x = -log(u)/lam
    return x

# Task 2

lam = 2
x = sampleExpo(lam,1000)

pylab.figure()
pylab.hist(x, normed=True)
xi = linspace(0,5,100)
pylab.plot(xi, lam*exp(-lam*xi))
pylab.show()

# Task 3

def sampleGammaInt(k, lam):
    return sum(sampleExpo(lam, k))

# Task 4

def gammaPDF(x, alpha, lam):
    return stats.gamma(alpha).pdf(x*lam)*lam

alpha = 5.7
k = floor(alpha)
lam = 2

M = gammaPDF(alpha-k,alpha,lam)/gammaPDF(alpha-k,k,lam-1)

xi = linspace(0,10,100)
pylab.figure()
pylab.plot(xi, gammaPDF(xi,alpha,lam),'b-')
pylab.plot(xi, M*gammaPDF(xi,k,lam-1),'r:')
pylab.show()

# Task 5

def sampleGamma(alpha, lam):
    k = floor(alpha)
    M = gammaPDF(alpha-k,alpha,lam)/gammaPDF(alpha-k,k,lam-1)
    while True:
        x = sampleGammaInt(k,lam-1)
        prob_of_accept = gammaPDF(x,alpha,lam)/(M*gammaPDF(x,k,lam-1))
        if random.random_sample(1)<prob_of_accept:
            return x

alpha = 5.7
lam = 2
x = [sampleGamma(alpha,lam) for i in xrange(1000)]

# Task 6

pylab.figure()
pylab.hist(x, normed=True)
xi = linspace(0,10,100)
pylab.plot(xi, gammaPDF(xi, alpha,lam))
pylab.show()