from numpy import *

import pdb
from numpy import random
#we COULD write these from scratch, but...
from PyMC import factln,fpoisson,fnormal
from PyMC import plot,show

#random number generators
rpois=random.poisson
rnorm=random.normal

from numpy.random import random

#some useful functions

def logit(p):
    x=log(p/(1.-p))
    return x

def expit(x):
    p=1./(1.+exp(-x))
    return p

#our posterior function

def posterior(p,x,n):
   # pdb.set_trace()
    like=x*log(p)+(n-x)*log(1.-p)
    return like
    
#metropolis-hastings    

def mh(start,iter,burn,x,n,tune):
    samples=0
    accept=0
    p=start
  
    p_=[]
    sd=tune
   
    while samples<iter:
        x_prop=rnorm(logit(p),sd)
       # pdb.set_trace()
        p_prop=expit(x_prop)
        
        J_prop=fnormal(x_prop,expit(p),1./sd**2)
        J=fnormal(expit(p),x_prop,1./sd**2)
        
        f_prop=posterior(p_prop,x,n)
        f=posterior(p,x,n)
        r= exp((f_prop-J_prop)-(f-J))
        if random()<r:
            p_.append(p_prop)
            p=p_prop
            accept+=1
        else:
            p_.append(p)
        
        samples+=1

    acceptance=float(accept)/float(samples)

    return p_,p_[burn:samples],acceptance

n=100
x=68

iter=10000
burn=5000

Full_trace,X,acceptance=mh(.5,iter,burn,x,n,.6)


iter_=arange(0,len(X))
plot(iter_,X)



print 'acceptance rate',acceptance
print 'variance', var(X)
samples=len(X)
X.sort()

lower=int(samples*.025)
upper=int(samples*.975)


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

print 'mean',mean(X)

show()