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

x=(0,2,4,6,8, 10, 12, 14, 16, 18,1,  3,  5,  7,  9, 11, 13, 15, 17, 19)

y=(1.4943758604867399,5.0499971697255219,9.604759217795511,13.086903746351343,18.539641921134745,23.180373381613105,27.497263381109164,32.946067205479302,38.000552082346204,42.551595817815411)


class m1(MetropolisHastings):

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

        
        self.parameter('b',init_val=[0.,2.0])
        self.node('y_miss',[10,])
        self.parameter('tau',init_val=1.0)
       
    def model(self):      
        
        self.normal_prior(self.b,0.,0.0001)
        sd=1./sqrt(self.tau)
        
        self.half_normal_prior([sd],[0.0001])

        


        
        mu=[self.b[0]+self.b[1]*z for z in x[0:10]]

        self.normal_like(y,mu,[self.tau]*len(y))

            
            
        for i in range(10):
            self.y_miss[i]=self.b[0]+self.b[1]*x[i+10]
        
            
        
        
        
         
        
  
def run(iter=20000,burn=10000,thin=10,model='m1',plot=False,conv=False,auto=False):

    #the instance of the sampler is stored in the dictionary
    models={'m1':m1()}
    #now we invoke sampler reps times, each time simulating the data with specified vectors n, p

    #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)

    sampler.export(results,'regression1')

    #replaces data object with new simulated data for next iteration


if __name__=="__main__":
    run()