# -*- coding: utf-8 -*-
## a rejection sampler to simulate values from the 
## beta density

from scipy import random
from scipy import stats
from scipy.special import gamma
from numpy import *
import pylab

## Madeleine's code for the Beta Density
def BetaDensity( x, a, b ):
	dens = gamma( a + b )/gamma( a )/gamma( b ) * x**(a-1) * (1-x)**(b-1)
	return dens


## specify the parameters of the Beta Density we wish to simulate
a     = 2.7
b     = 6.3

## estimate the density at the maximum of our chosen distribution
maxd  = BetaDensity(1.7/(1.7+5.3),2.7,6.3)

## simulate 1000 candidate values
x     = random.uniform( 0,1,1000 )
## simulate 1000 uniforms
u     = random.uniform( 0,1,1000 )

## calculate f(x)/Mg(x)
test  = BetaDensity(x, a, b)/maxd*1
## note that the density of the uniform is 1 everywhere


## create an array of accepted values
out   = x[u<test]


## Plot a histogram of the results
pylab.figure(1)
pylab.hist(out,bins=10,range=[0,1],normed=1)
##pylab.show()



## Use Madeleine's code to plot the beta density
pivec = linspace( 0, 1, 100 )
pivec = pivec[2:99]

bdvec = BetaDensity( pivec, a, b )
pylab.figure(1)
pylab.plot(pivec,bdvec,color=(0,1,0))
pylab.title('Beta density')
pylab.xlabel("x")
pylab.ylabel("Density")
pylab.show()