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




def simdata(n,p):

    #now we'll simulate the data

    success=[binomial(w,z) for w,z in zip(n,p)]
    trials=n
    reps=len(n)
    return trials,success,reps




class m1(MetropolisHastings):

    def __init__(self,n,p):     
        'Initialization of model parameters and nodes' 
        MetropolisHastings.__init__(self)

        self.data=simdata(n,p)
        self.parameter('p',init_val=[0.5]*self.data[2])
       
    def model(self):      
        n,x,reps=self.data
        self.uniform_prior(self.p,0,1)        
        self.binomial_like(x,n,self.p)
        
        
        
         
        
  
def run(iter=20000,burn=10000,thin=10,model='m1',plot=False,conv=False,auto=False,p=[0.5]*5,n=[100]*5,reps=2):

    #the instance of the sampler is stored in the dictionary
    models={'m1':m1(n,p)}
    #now we invoke sampler reps times, each time simulating the data with specified vectors n, p
    for i in range(reps):
        #first time through model self initializes with simulated data
        sampler=models[model]
        sampler.sample(iterations=iter,burn=burn,thin=thin,plot=plot)

        results=sampler.summary(burn=burn,thin=thin)
        if conv:
            sampler.convergence(burn=burn,thin=thin)
        if auto:
            sampler.autocorrelation(burn=burn,thin=thin)
        name='stats'+str(i)
        sampler.export(results,name)

        #replaces data object with new simulated data for next iteration

        sampler.data=simdata(n,p)

if __name__=="__main__":
    run()