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

x=(0,2,4,6,8, 10, 12, 14, 16, 18)


y=(1.4943758604867399,5.0499971697255219,9.604759217795511, 13.086903746351343, 18.539641921134745, 23.180373381613105, 27.497263381109164, 32.946067205479302, 38.000552082346204, 42.551595817815411,2.3818651120545664, 8.9120638434691877, 10.918265656363822, 17.759741682053594, 22.197026902672608, 26.402869557805943, 30.584302148927126, 34.488676852266366, 41.25177717156442, 43.465687300059443)

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('x_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]

        self.normal_like(y[0:10],mu,[self.tau]*len(x))

    

            
            
        for i in range(10):
            
            self.x_miss[i]=(y[i+10]-self.b[0])/self.b[1]
        
            
        
        
        
         
        
  
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,'regression2')

    #replaces data object with new simulated data for next iteration


if __name__=="__main__":
    run()