from numpy import *
from numpy.random import random
import pdb


#rejection sampling example
#uniform envelop, exponential targe

mu=0.56

samples_needed=10000
tries=0
samples=0
range_=[0,10]
X=[]
while samples<samples_needed:
    x=random()*range_[1]+range_[0]   #takes uniform random value in range
    g=mu   #forces uniform to be >= f(x)

    f=mu*exp(-mu*x)  #calculates f(x)
    alpha=f/g       #importance ratio
    tries+=1
    if random()<alpha:
        X.append(x)
        samples+=1

print 'number of samples', samples
print 'acceptance rate',float(samples)/float(tries)
print 'mean', mean(X)
print 'variance', var(X)
X.sort()
lower=int(samples*.025)
upper=int(samples*.975)


print 'lower',X[lower]
print 'upper',X[upper]
print 'median',median(X)