# -*- coding: utf-8 -*-
"""
Created on Thu Feb 16 11:57:21 2017

@author: camacho
"""

import random,math
 
def gibbs(N=1000,thin=500):
    x=0
    y=0
    print "Iter  x  y"
    xresul=[];yresult= []
    for i in range(N):
        for j in range(thin):
            x=random.gammavariate(3,1.0/(y*y+4))
            y=random.gauss(1.0/(x+1),1.0/math.sqrt(x+1))
            xresul.append(x);yresult.append(y)
        #print i,x,y
    return [xresul,yresult]
    
a=gibbs()